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I. 


INTRODUCTION 


i 


Structural  failure  of  reinforced  concrete  structures, 
as  well  as  other  kinds  of  structures  must  be  defined  in 
terms  of  the  degree  of  damage.  In  this  study  failure  is 
defined  as  the  actual  fracture  or  separation  of  a reinforc- 
ing element  in  the  slab  or  beam.  This  means  that  before 
failure  occurs  considerable  concrete  cracking  will  have 
occurred  and  the  reinforcing  element  will  have  yielded  and 
reached  the  ultimate  stress  or  strain.  Since  the  tensile 
strength  of  concrete  is  negligible  in  comparison  to  the 
tensile  stx’ength  of  the  reinforcing  element  the  failure 
criterion  is  based  solely  on  the  ultimate  stress  or  elonga- 
tion of  the  reinforcing  elements . 

In  this  study  reinforced  concrete  beams  and  plates  sub- 
jected to  transverse  blast  or  impulsive  loads  are  considered 
for  analysis.  Experimental  investigations’*'  cf  two  edges 
fixed  reinforced  concrete  slabs  subjected  to  small  under- 
ground near  field  explosions  showed  "plastic  hinge"  mecha- 
nisms at  failure.  Typical  failures  of  these  type  are  shown 
in  Figure  1.  Failure  of  this  type  suggest  the  same  type 
failure  mechanism  as  that  observed  in  metal  beams  when  sub- 
jected to  impulsive  loads  of  magnitude  greater  than  the 
static  collapse  load  of  the  member.  For  plates  and  beams 


_ • . . • > w, 


a)  36"x36"x4"  two  edges  fixed 
8.0  lbs  @ 3 ft. 


cal  failures  o 


with  transverse  or  normal  loadings  where  bending  is  pre- 
dominant an  idealized  rigid  perfectly  plastic  constitutive 
relation  can  be  used  to  estimate  the  true  moment  carrying 
capacity  of  the  structure.  This  assumption  appears  justi- 
fied for  the  ductile  steel  reinforcing  elements  whose  stress 
strain  curves  show  a rather  small  elastic  portion  when  com- 
pared to  the  overall  stress-strain  curve  for  the  material. 
The  "fully  plastic"  or  "ultimate  moment"  hinge  as  applied 
to  statically  loaded  reinforced  concrete  structures  is 
summarized  by  Szilard^.  The  "plastic  hinge"  or  "fully 
plastic  moment"  is  covered  in  detail  by  Timoshenko  and  Gere 
Abrahamson,  et  al,1*  used  both  stationary  and  moving  "plastic 
hinges"  to  predict  response  of  thin  metal  beams  and  plates 
when  subjected  to  spatially  uniform  loads  of  several  differ- 
ent time  variations.  Jones"*  used  this  simple  method  to 
successfully  predict  the  large  inelastic  behavior  of  thin 
metal  beams  when  subjected  to  impulsive  loadings  and  found 
very  good  correlation  with  the  experimental  results  of 

6 7 Q 

Menkes  and  Opat  . Both  Florence'  and  the  authors  have 
applied  this  method  to  reinforced  concrete  beams  and  rea- 
sonable agreement  was  found  when  analytical  results  of 

Reference  8 were  compared  to  the  experimental  results  of 

• fi 

1.  Based  on  the  results  obtained  by  the  authors 

3 . 


Reference 


r. 


for  beans  with  linear  spatial  leadings  the  study  was  ex- 
tended to  cover  both  linear  and  nonlinear  impulsive  spatial 
loadings  for  fixed  and  simply  supported  reinforced  concrete 
beams  and  plates. 

II.  ANALYSIS 

2.1  Assumed  Response  Mechanisms 

The  dynamic  response  of  reinforced  concrete  plates 
and  beams  may  well  be  broken  down  into  four  failure  mech- 

[ 

anisms  as  described  below.  In  each  particular  case  the 
failure  mechanism  exhibited  by  the  structure  is  dependent 
on  the  initial  magnitude  of  the  assumed  loading  or  pressure 
function.  It  is  assumed  that  no  response  will  occur  unless 
the  initial  value  of  the  loading  pulse  has  exceeded  the  static 
collapse  load  of  the  structure  based  on  a rigid  perfectly 
plastic  constitutive  relation.  The  static  collapse  load  is 
that  load  which  will  cause  initial  response  of  the  plate  or 
beam  in  the  mechanism  1 mode  described  below. 

Response  under1  dynamic  leadings  may  be  characterized 
by  the  following  four  mechanisms. 

Mechanism  1.  This  response,  as  shown  in  Figure  2, 
occurs  when  thm  initial  magnitude  of  the  .loading 
pulse  is  slightly  larger  than  the  static  collapse 
loac  tor  the  member. 

*+ . 


& 


!)!C' 
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Figure  2.  Failure  by  a mechanism  1 mode. 


Mechanism  2.  As  the  initial  magnitude  of  the  loading 
pulse  increases  a plastic  hinge  at  interior  points 
'other  than  at  the  center  of  the  member  will  form  and 
move  toward  the  center  of  the  member.  This  mechanism 
type  for  a beam  is  shown  in  Figure  3. 


M-V, 

O I 
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Figure  3.  Mechanism  2 type  beam  response. 


Mechanism  3.  As  the  initial  magnitude  continues  to 
increase  the  initial  hinge  location  x^  moves  closer 
to  the  end  of  the  beam  and  a discontinuity  will  exist 
in  the  analysis  as  xj,  = 3.  This  means  failure  by 
complete  shear  at  the  edges  will  occur  before  any 
hinge  motion  can  occur.  This  type  failure  for  metal 
beams  has  been  discussed  by  Jones^,  however  a general 
discussion  of  failure  of  reinforced  concrete  beams 
by  this  mechanism  is  not  included  in  this  report. 

Mechanism  4.  For  very  local  blast  loads  on  long  beams 
or  large  plates,  failure  of  the  concrete  begins  before 
any  appreciable  overall  response  can  occur  and  the 
concrete  and  reinforcing  elements  fail  locally.  This 
type  mechanism  is  very  similar  to  that  of  mechanism  3 
above  and  discussion  of  this  mechanism  is  also  not  in- 
cluded in  this  report. 

For  "plastic  hinge"  response  in  a static  case,  the 
assumption  is  made  that  when  the  moment  at  a point  in  the 
beam  reaches  sufficient  magnitude  to  cause  yielding  of  the 
reinforcing  elements  then  rotation  will  increase  unbounded 
to  failure  unless  the  load  is  removed  in  time  to  prevent 
total  failure.  This  type  failure  mechanism  may  be  described 
in  terms  of  a uniform  static  load  on  a beam  as  shown  in 
Figure  4. 

6 . 
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P (force/unit  length) 


MOMENT 


SHEAR 


Figure  4.  Load,  shear  and  moment  diagram  fo 
uniformly  loaded  (static)  fixed  ended  beam. 


The  maximum  positive  bending  moment  occurs  at  the  midspan 


and  the  maximum  negative  bending  moment  occurs  at  the  ends 


of  the  beam.  Rotation  at  the  "hinge"  begins  when  this 


ufficient  to  cause  yielding  in  the  reinforcing 


moment  is 


may  be  determined  from 


element.  This 


moment 


ection  of  the  beam  as  shown  in  Figure  5,  and  is 


reinforcing  element 
in  tension 


Figure  5.  Beam  cross  section  for  calculating  a 
plastic  moment  MQ. 


based  on  the  assumption  that  the  tensile  strength  of  the 


e 6.  This  moment  is 


The  hinge  moment  per  unit  length  is  defined  as 


b are  shown  in  F 


Figure  6.  Assumed  stress  strain  curve  of  the 
reinforcing  elements. 

stress  of  the  reinforcing  element,  oc  , is  the  compressive 

strength  of  the  concrete,  and  q is  the  reinforcing  ratio, 

i.e.,  the  area  of  reinforcing  elements  in  tension  to  the  < 

total  cross  sectional  area  of  the  beam. 

When  the  moment  at  the  midspan  reaches  M , localized 
rotation  occurs  at  the  ends  and  midspan  and  the  resulting 
response  is  as  shown  in  Figure  2.  For  a static  case  rota- 
tion will  occur  until  fracture  of  the  reinforcing  element 
occurs.  Rotation  for  failure  of  the  reinforcing  elements 
will  be  discussed  later.  The  remainder  of  this  report  is 
restricted  to  discussion  of  mechanisms  1 and  2 dynamic 
response  in  beams  and  plates  when  subjected  to  spatially  and 
time  varying  loading  functions.  One  very  important  and 
underlying  assumption  used  in  this  analysis  is  that  the 


shear  at  all  interior  hinges  of  the  beams  and  plates  is 


* 

I 

i 


s 


i 
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zero.  This  is  based  on  the  observations  of  the  shear  diagram 
of  Figure  4 and  is  discussed  by  Abrahamson4  for  both  the 

mechanism  1 and  2 modes. 


2.2  Assumed  Loading  Functions 

The  loading  functions  used  to  simulate  the  blast  or 
underground  explosion  is  simplified  by  breaking  it  down 
into  two  parts  with  separable  time  and  spatial  variables. 
The  general  expression  for  the  loading  is  given  as 

P(x,y,t)  = f(t)P(x,y)  (3) 

where  x,y  are  the  spatial  coordinates  and  t is  the  time. 

.Two  types  of  time  variations  were  chosen  for  the 
analysis.  The  fundamental  square  wave  or  impulse  of  the 
form 


f(t)  = 1 , 0<t<x 

f(t)  = 0 , t>x 

and  a decay  type 

f(t)  = (1-t/i  )exp(-at/x  ) CKt<_i 
f(t)  = 0 , t > T 


(4) 


(5) 


where  x is  the  pulse  duration  in  time  and  a is  the  time 
decay  constant.  These  two  forms  are  shown  schematically 
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a)  square 
Figure  7. 


wave  b) general  decay 

Schematic  of  time  function  f(t). 


in  Figure  7.  The  triangular  pulse  may  be  obtained  from 
Equation  (5)  when  a = 0.  No  difficulty  arises  when  a -*■  0 
as  long  as  the  time  integrals  are  accomplished  numerically, 
however  if  a general  closed  form  solution  is  required  then 
the  triangular  or  linear  decay  time  function  must  be  handled 
as  a special  case.  This  problem  arises  due  to  the  1/a  appear- 
ing in  the  result  of  integrating  Equation  (5)  with  respect 
to  time. 

In  this  study  the  spatial  integration  was  handled  in 
closed  form  and  the  spatial  loading  function  P(x,y)  was 
assumed  to  be  of  the  symmetric  form  of  Figure  8 and  given 
as 


P(x) 


P(y ) 


0<x<a 


0<y<b 
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(6) 


. i 


b ) bean: 


for  beams  where  PF  is  a uniform  load  acting  ove 


the  w’hole 


non 


showing  the  coordinates  (x,y)  and 


a,b,  of  both  plate  and  beam  spatial  loadings  are 


dimension 


cases  of  Equations  (5)  and  (6)  where  3=0 


integration  was  handled  in  closed  form  and  a separate  case 
for  3=0  was  included  in  the  analysis.  In  Equation  (5)  the 
same  decay  constant  3 was  used  for  both  the  x and  y variations. 
This  would  be  true  only  for  a central  symmetric  load  on  a 
square  plate.  However,  due  to  uncertainities  in  both  the 
mechanics  model  and  loading  functions  the  assumption  appears 
reasonable  for  aspect  ratios  of  at  least  two  or  less.  Both 
the  spatial  and  time  decay  terms  a and  3 are  obtained  from 
observations  of  either  experimental  or  analytically  determined 
pressure  time  histories  of  under  ground  blasts. 

2.3  Equations  of  Motion 
2.3.1  Introduction 

The  equations  of  motion  for  both  beams  and  plates  are 
based  on  a conservation  of  energy.  The  equations  as  presented 
here  were  derived  from  force  balance  equations  and  verified 
by  comparing  with  independent  calculations  on  individual 
energy  terms.  Beams  whose  response  starts  with  interior 
hinges  not  at  the  center  (mechanism  2)  are  assumed  to  respond 
in  this  manner  until  the  moving  hinges  reach  the  beam  center 
and  respond  in  a mechanism  1 mode  until  the  rotation  stops. 

For  plates  the  location  for  the  mechanism  1 response  is 
dependent  on  both  the  load  shape  and  plate  aspect  ratio  and 
may  be  determined  by  a method  described  by  Szilard"*.  For 


plates  starting  in  a mechanism  2 response,  it  is  assumed  that 
the  final  hinge  locations  will  be  that  determined  as  if  the 


I 


plate  assumes  a mechanism  1 mode.  This  means  that  for  each 
plate  a final  hinge  location  must  be  calculated  for  each 
loading.  The  method  for  doing  this  is  discussed  along  with 
the  derivation  of  the  equations  of  motion  for  the  plate. 

2.3.2  Beams 

Since  the  beam  response  is  either  mechanism  1 or 
mechanism  2 with  interior  hinges  moving  toward  the  center, 
the  derivation  will  be  based  on  a mechanism  2 mode  and  the 
mechanism  1 response  will  be  a special  case  when  the  hinges 
start  at  the  midspan  of  the  beam.  The  beam  loading  and  response 
is  assumed  to  be  completely  symmetrical  therefore  only  half  the 
beam  need  be  considered.  Using  the  response  shape  of  Figure 
9 the  derivation  of  the  beam  equations  of  motion  is  as  follows. 
For  the  exterior  portion  of  length  x^  the  general  equation  is 


10 


P(x,t)xdS-n 


wxdS-FMQ 


(8) 


where  I is  the  moment  of  inertia  and  0 the  angular  accelera- 
tion about  the  axis  at  0,  w is  the  weight  per  unit  area  of 
surface  area,  S is  the  surface  area  over  which  the  pressure 
is  acting,  n is  a position  indicator  with  n = 0 for'  a vertical 
wall,  n=-l  for  bl  ,t  above  beam  and  n~  + l for  blast  below  beam, 
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Figure  9.  Schematic  for  beam  response  of 
mechanism  2. 


and  F is  an  end  fixity  term  with  F=2  for  fixed  ended  beams 
and  F=1  for  simply  supported  beams.  Using  Equations  (2) 
and  (7)  with  m the  mass  per  unit  of  surface  area,  w the 
weight  per  unit  of  surface  area  and  dS=b  dx  Equation  (8) 
becomes 


Performing  the  indicated  integrals  and  with  X=x^/a  Equation 


(9)  becomes 


-----  t-)  ■ (pJ^Pc  [(  B2 X2  - 2 BX+  2 ) expB (X-l ) 
3max3  L 

-2exp(-B)]l 


3nw  31 Mu 
2maX  ma 3 X^ 


For  the  inner  portion  of  length  (a-xp)  the  general  equation 


6 bm  dx  = f(t)  bLPv+Pp  ^lexpB  1 1 ] dx  -n  bwdx  (11) 

Lj  ^ \a)  U J 


Performing  the  indicated  integrations  Equation  (11)  becomes 


f(t)  r 

\ n_..o2  . 


<5  = j pEe  +i3x[e"1_(ex~1)exP3(x“1)]j 


In  order  to  maintain  continuity  between  the  two  elements 


9xp  = 0aX  = 6 


8xv  = 0aX  = 6 


The  assumption  is  made  that  Equat ions  (10),  (12)  and  (13) 
satisfy  all  the  conditions  required  for  motion  by  plastic 


hinges  and  in  order  to  determine  the  initial  response  the 
initial  conditions 


l 


0aX 
f (t) 


(14) 


are  imposed  on  equations  (10)  and  (12).  The  resulting 
equation  becomes 


a2  B 3 P-pY*  - nw  3 3 a * Y * - 6FM,  3; 
*-•  u 


3 _ 2 v 2 


(lb) 


+ 2a2Pc-Ci-[l-3+(3Y-l)exp3(Y-l)]+3[(32Y2-23Y+2)exp3(Y-l) 

L 1 

-2 exp (-3)  JJ 


= 0 . 


Where  Y,  the  initial  value  of  X(xft0/a)  determines  the 
mechanism  for  initial  response.  Only  solutions  for  0<Y<1 
are  allowed  and  solutions  where  Y>1  are  set  Y=l.  The  solution 
Y=1  means  a mechanism  1 response  and  for  this  case  x^  is  set 
equal  to  a in  Equation  (10)  and  the  response  is  determined 
from  this  equation  only  and  Equation  (1?)  is  disregarded. 

For  the  case  of  a linear  spatial  variation  for  the  beam,  i.e., 
if  3=0  in  Equation  (7),  then  the  equation  of  motion  must  be 
derived  for  this  case  and  the  0,5  equations  similar  to  Equa- 
tions (10)  and  (12)  are  found  to  be 


* 


e = 


f(t) 


2 maX 


y(3Pj,+  2XPc)  - 


3nw 

2maX 


3FM, 


(16a) 


ma ' 
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6 = ^^-fpE+— (X+l)]  . (16b) 

rn  v L ? J m 

Applying  the  same  initial  conditions  as  in  Equation  (14)  and 
the  same  procedure  for  finding  the  initial  response  of  Equa- 
tions (10)  and  (12)  then  the  beam  response  for  3 = 0 may  be 
determined. 

Assuming  conservation  of  energy  the  equations  of  motion 
can  be  double  checked  using  the  relation 

• • • • 

KE  = WF  - WP  - PE  (17) 

where  the  (')  means  time  derivative,  KE  is  kinetic  energy, 

WF  work  done  by  applied  pressure  or  force,  WP  work  done  by 
plastic  hinge  and  PE  is  potential  energy.  In  that  it  is 
desirous  to  monitor  these  values  for  given  response  calcula- 
tions, terms  for  the  right  hand  side  of  Equation  (17)  were 
derived.  For  the  nonlinear  spatial  function  (8^0)  the  terms 
are 


WF  = 


WP  = 


2ba2  &f ( t ) f y 

— |Pe33X(1-|)  + PcC(2-6X)exp3(X-l)-2exp(-3) 

+3X(3-1)]J 


(18) 


2bFMuG 


PE  = 2a2bnwX(l-|)9 


where  w is  the  weight  per  unit  area.  The  work  and  energy 
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terms  for  the  linear  spatial  loading  (3=0)  are  given  as 


WF  = 2a2bf  ( t ) 0 [PEX(1-|]  + ^X  [p£j] 


WP  = 2bFMu0 
PE  = 2a2bnwX 


(19) 


The  kinetic 
integrating 

KE  = 


energy  term  is  calculated  by  numerically 
Equation  (17)  to  yield 

(WF-WP-PE)dt 


( 20) 


2.3.3  Plates 

For  plates  the  hinge  positions  for  static  loading  are 
dependent  on  plate  size  and  loading  and  must  be  calculated 
for  each  type  of  load.  This  method  used  here  is  an  exte.n- 

• • 9 

sion  of  the  method  used  by  Szilard4  for  uniformly  loaded 
plates.  For  square  plates  and  uniform  loadings  the  hinges 
of  the  plate  are  assumed  to  run  from  each  corner  to  the  plate 
center  as  shown  in  Figure  10a.  For  rectangular  plates  and 
non  uniform  loadings  the  hinges  are  assumed  to  form  as  shown 
in  Figure  10b.  The  interior  hinge  A-B  of  Figure  10b  is  always 
assumed  to  run  parallel  to  the  long  side  of  the  plate. 
Schematics  of  the  deformed  shape  are  shown  along  with  the. 

19. 


quarter  of  the  plate  shown  in  Figure  11.  Justification  of  this 


yinmetric  response  shape  and  loading 


(hinge  line) 


Figure  11.  Schematic  used  for  calculating  final 
hinge  locations. 


model 


a given  unit  center  displacement  yields  an  equation  in  terms 


of  the  plate  dimensions,  loading  parameters,  and  the  inter 


zb.  For  the  derivation,  0 is  assumed  to 


be  the  rotation  along  the  line  OD  and  the  assumption  of  con 


tant  displacement  along  BC  the  rotation  about  line  OE  is 


placement  along  the  line  EB  is  then  6 = y/zb  and  the 


displacement  along  the  line  lC  is  6 = x/a.  The  xternal 
work  done  by  the  pressure  load  of  Equation  (6)  over  the  area 
S is  in  general 


WeXf  = P6  dS 


and  upon  integration  over  the  quarter  of  the  plate  and  mul- 
tiplying by  four  gives 


. = 4ab 
ext 


h(f9 


where 


f(3»z)  = ( 3-0z)exp£  ( z-1 ) - ( 2gz2 + 9z2  + 2flz+ 3 ) exp  ( - 3 ) 
+262-73z2+9z2  . 


The  internal  work  is  the  work  done  by  the  plastic  hinge 
moving  through  some  local  angle  at  the  hinge  lines.  The 
work  done  by  the  hinge  lines  that  move  diagonally  are  based 
on  the  projected  component  on  the  coordinate  axes.  For  in- 
stance the  line  OB  of  Figure  11  does  work  an  amount  of 
(Mua)(6a/zb)  about  an  axis  parallel  to  ox  and  an  amount 
of  (Muzb)('i)  about  an  axis  parallel  to  the  oy  axis.  Based 
on  this  assumption  and  depending  on  the  edge  conditions 
the  internal  work  may  be  expressed  as 


Wint  = 


"-GH) 


where  0a  = 1 for  the  assumed  displacement  A = 1.  Equating 


Wext  to  W^nt  and  rearranging  terms  yields 


Ppb2 


(z+AR2z2)34 


f (3,z)+K04z 


4_2 a z 


P g 

» K = p—  . 

* C 


If  the  left  hand  side  of  Equation  (25)  is  plotted  against 


z a minimum  value  is  realized.  It  was  found  analytically 


that  only  one  minimum  existed  for  0<z<l  therefore  the  value 


of  z for  a minimum  can  be  determined  by  differentiating 


Equation  (25)  with  respect  to  z and  setting  the  right  side 


to  zero.  The  results  of  this  for  the  nonlinear  spatial 


load  is 


AR2z[p£04z2  f|)+PCC  (6-43z  + g2z2)exp3(x-l)-(6  + 23z)exp(-3)]j 
+ PE34  z2^y-^j  +PC[ (3-33z+32z2) exp 3 (z-l)  + (25z2  + 9z2-3)exp(-3)] 


-232z2+73z2-9z2  = 0 , 


and  for  the  linear  spatial  load 


P£ ( 4AR  2 z+  8z- 12 ) + PC ( 2AR2  z3  + 3z2-5)  = 0 


The  solution  of  Equations  (26)  or  (27)  for  z determines  the 
hinge  locations  for  a mechanism  1 response  and  final  locations 
for  the  mechanism  2 response  in  plates. 

The  schematic  of  a mechanism  1 plate  response  would  be 
identical  to  that  of  Figures  10  and  11  and  the  schematic  of 
a mechanism  2 plate  response  is  shown  in  Figure  12. 


a)  3D  schematic  b)  top  view 


Figure  12.  Mechanism  2 plate  response. 

As  in  the  case  for  the  beams,  only  the  derivation  of 


the  mechanism  2 response  will  be  given  and  when  x,,  = a only 
the  equation  of  motion  for  6 will  be  used.  For  this  deriva- 
tion many  of  the  same  terms  used  for  the  beam  will  also  be 


used  here.  Most  of  these  terra 


are  shown  in  Figure  13  and 


only  a quarter  of  the  plate  is  shown  due  to  symmetry  of 


loading  and  response.  The  assumption  for  mechanism  2 is 


bx  (pressure 
a line) 


Schematic  for  derivation  of  equations 
plate. 


that  the  inner  portion  of  the  plate  is  flat  and  has  a displace- 
ment 6,  same  as  that  of  the  plate  center  point.  This  in- 
sures continuity  of  the  inner  hinge  line  and  as  shown  in 
Figure  13a  gives  different  rotations  for  the  plate  edges 
when  z < 1.0.  Also  for  the  case  of  z < 1.0  the  distances 
between  the  moving  hinges  and  their  respective  edges  are 
not  equal.  However,  for  the  square  plate  with  z = 1.0 
both  edge  rotations  will  be  equal  and  moving  hinge  locations 
will  be  the  same. 

Using  Figure  14  the  integrals  for  equation  of  motion 
were  broken  into  five  separate  areas  (T)  to  (5)  . For  the 
outside  areas  the  equation  of  motion  is 


Figure  14.  ' use  Lr  integrals  of  equation 

of  motion. 


?•= . 


using  AK  = b/a,  X = x^/a,  evaluating  the  integrals  and  re- 
arranging terms  Equation  (28)  becomes 


m0a4AR 


<X3  |c  z 2AR+1 ) fi-2T)  + ( 1- z 
L v3  hJ  4. 


_ f (t)a3?HT 
z234 


-nwa 3 ARX 2 ARz2 


^ [pe34z2X2[~ARz2  [l-|j+i(i-z)x+l] 

+PcC-( 2AR3z2+20z2+3z2+5ARz-3)exp(-B) 

+(g3z2X2-83z2X3+432z2X2) exp 3 ( X-l ) 

+(-232z2X-93z2X+23z2+9z2 )exp3(X-l) 

+ AR  (-33z4X3  + 33z,  '*X2  + 332z3X2 

-2  8%;3X-63z2X+23z2  + 6z)exp3(  zX-1 ) 

+ ( ~32  z2X2  + 3gzX-3 )exp3 ( zX-1 ) ]J 

ARz  2 [---!+  ( 1-z )— +i~— 1 
L U 3J  3 2 3J 


FMu(l+AR)a 


For  the  inner  flat  portion  the  integrals  are  broken 
into  2 parts,  i.e.,  areas  (V  and  (T)  . In  this  case  the 
interia  term  is  much  simpler  and  is  included  under  the  same 
integral  as  the  potential  energy  term. 


(m<5  + nw)  (a-x^)  |^b-— XjJ 


f (t) 


bx 
■ a 

(pe*pc 


xh 


expB 


dy  dx 


(30) 


Evaluating  the  integrals  and  rearranging  terras  Equation  (30) 
becomes 


m(l-X) (l-zX)6  = f (t)JpE(l-zX+zX2-X) 

PC 

+— [ 23-4+(  3 + 32zX2-32zX-3)exp3(zX-l) 
e3  (31) 

+ ( 4 + 3-3  3X-3  ZX+  3 2X2 )exp3 (X-l ) ]J 


-nw( 1-X ) ( 1-zX) 


Using  the  same  approach  as  in  the  case  of  the  beams 
the  kinetic  energy,  plastic  work,  potential  energy  and 
pressure  work  of  the  plate  with  nonlinear  spatial  loads 
(0^0)  may  be  calculated  using  the  following  equations. 


:>  , 3 


_ 4 0 f ( t ) ARa 


2 ok 


Z 3 


h-'-t'-f'f  t) 


+Pc[23zzzX-43z2X) 


+(-32z2X+B2z2X2-33zX+23z+3)exp3(zX-l) 


+(-53z2X+32z2X2-32z2X+23z2+9zz)exp3(X-l) 


(23z2+23z+9z2+3)exp(-3)]^ 


(32) 


PE  = 4nwARa  0 X 


3fllv  X2  zX2. zX3 


0 [x 


2 2 3 


-) 


WP  = 4 F0M 

For  the  linear  spatial  loading  (3=0)  the  sane  procedure 
as  above  was  followed  in  deriving  the  equations  of  notion 
for  the  plate.  Again,  the  reason  for  the  separate  deriva- 
tion for  3=0  is  that  in  the  nonlinear  spatial  loading  equa- 
tions such  as  Equations  (29)  and  (31)  the  3 term  appears  in 
the  denominator . In  these  equations  L ^Hospital  rule  could 
be  used  as  O^O , however  it  appeared  to  be  much  easier  and 
straight  forward  to  simply  derive  the  equations  from  basic 
assumptions.  With  these  comments  in  mind  and  using  the 
same  general  procedure  leading  to  Equations  (29)  and  (31) 

30  . 


- , _ , 


■ 


the  equations  of  motion  for  a plate  with 


0 are 


3 f / _ 2 


m6a  ARX 3 [ ( z2AR+l)  (j~£)+  ( 1-z )-] 


f(t)a3ARX2^PE  J(z2AR+1)[j-|J  + (1-X)|"| 


+ PC  [z3XAR  f— — 'l  + — 

L L I3  y 4 


(t-t-)]  } 


-nwa3 ARX2 


[AR<H>(1-2tH] 


-FMu(l+AR)a  , 


m(l-zX)  0=  f(t){PE(l-zX)+Pc  j(X+l)-^-^z2+jj]  ]• 


-nw(l-zX)  , 


with  continuity  assured  by  0Xa  = 6 . The  corresponding  equa- 


tions for  WP  and  PE  are  the  same  as  Equation  (32)  except  for 


WF  which  is 


WF  = 4f(t)9ARa3  Pr  X 


fp4x 


X2  zX2  zX 3' 


2 2 3 


+ Pp 

c U 6 6 6 


6 6 4 


-^)]  • 
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The  general  solution  procedure  here  is  the  same  as  for 
beams.  After  the  inner  hinge  position  z has  been  determined 
by  either  Equation  (26)  or  (27)  it  becomes  a constant  for 
the  equations  of  motion.  In  order  to  determine  the  initial 
positions  for  the  inner  hinges  the  initial  conditions 


6aX  = 6 


f(t)  = 


t = 0 


are  imposed  on  the  Equations  (29)  and  (31)  or  Equations  (33) 
and  (34).  Due  to  the  length  and  complexity  of  these  equa- 
tions they  were  omitted  but  included  in  the  numerical  calcula- 
tions to  be  discussed  later.  The  solution  resulting  from 
these  initial  conditions  determines  the  initial  mechanism  to 
be  used.  Again  assuming  Y is  the  initial  value  of  x^/eCx^/a), 
then  for  0<Y<1  a mechanism  2 response  is  initiated  and  the 
response  is  governed  by  both  the  6 and  6 equations  until 
the  inner  hinges  reach  the  final  positions  with  Xft=a.  Also 
for  Y=1  the  initial  rc sponse  is  a mechanism  1 type  and  only 
Equations  (29)  or  (33)  d be  applicable  for  a plate  problem. 

This  conclude.  t..e  :•  rivation  of  the  equations  of  motion 
necessary  for  de.  \ err  ining  the  dynamic  response  under  the 
assumed  pl  stic:  hinge  motion.  The  range  of  applicability  of 
these  equation  and  their  shortcomings  will  be  discussed 
later  under  a separate  section  but  before  that  a section 
ci  Lbing  a proposed  failure  criterion  will  be  given. 


I . 


2.4. 


Failure  Criterion 


Under  the  assumptions  used  in  the  previous  sections 
the  rotations  are  assumed  to  be  highly  localized  and  mathe- 
matically the  hinges  are  assumed  to  have  zero  width.  In  a 
real  or  true  structure  the  hinge  width  or  damaged  area  must 
be  non  zero  due  to  the  finite  thickness  of  the  structure. 
Intuitively,  the  thicker  the  structure  the  greater  the  hinge 
width  due  to  the  assumed  bending  taking  place.  For  an  assum- 
ed damaged  or  deformed  width,  at  a hinge  location,  the  total 
strain  may  be  estimated. 

The  tensile  strain  in  the  reinforcing  element  at  a 
hinge  position  may  be  broken  into  two  parts,  i.e.,  the  strain 
due  to  bending  and  the  strain  due  to  axial  elongation  in 
tension  e^.  The  total  strain  will  then  be  the  sum  of  these 
or 

e = + et  . (37) 

For  a given  assumed  damage  or  hinge  width  £ the  strain  terms 
may  be  approximated  using  the  diagram  of  Figure  15.  Using  a 
known  rotation  9 and  an  arc  length  £.  the  radius  of  curvature 
R at  the  hinge  becomes 

R = £ (38) 

a 

and  assuming  strain  due  to  bending  as  proportional  to  the 

33. 


Reinforcing  Elements 


Figure  15.  Schematic  of  plate  or  beam  thickness 
showing  rotation,  placement  of  reinforcing  elements 
and  assumed  hinge  length. 


distance  from  the  assumed  neutral  axis  then 


where  d is  the  distance  between  the  tensile  reinforcing 
element  and  outer  compressive  surface. 

For  a mechanism  1 mode  the  strain  due  to  axial  tension 
may  be  approximated  by  the  following.  Using  the  sketch  of 
Figure  16  the  increase  in  the  halfspan  length  A may  be  approxi 


Expanding  the  radical  in  a binomial  series 


Assuming  half  the  increase  A take 


uniform  over  the  deformed  length  £ then 


Approximating  the  rotation  as  6=6/a,  which  yields  6- Da,  and 
using  this  along  with  Equations  (41)  and  (42)  results  in 
the  strain 

et  = — (43) 

4H 

Using  Equations  (39)  and  (43)  the  total  strain  e becomes 


0a  fy  _cp| 

l (4  2 a J 


Examining  Equation  (44)  in  terms  of  the  ultimate  strain  of 
the  reinforcing  elements  and  solving  for  0U  required  to 
produce  ultimate  strain  gives 

0U  = /4£dcu 

a L/  "d  2 +1  -ij  . (45; 

For  a mechanism  2 response  the  same  argument  may  be  used 
and  the  resulting  equation  for  this  response  would  be 
identical  except  the  half  span  a would  be  replaced  with 
x^  and  the  Equation  (45)  for  a general  case  becomes 


_ d f /UxftEu  ^ 

6u  - +1  -5 


III.  DISCUSSION  AND  RLSULTS 


3.1  Introduction 

The  previous  sections  contain  the  necessary  equations 
of  motion  to  determine  the  dynamic  response  of  beams  and 
plates  for  spatially  and  time  varying  loading  functions  for 
the  assumed  plastic  hinge  motion.  The  equations  as  derived, 
with  adjustments  of  hinge  moment  term,  would  apply  to  thin 
metal  structures  as  well  as  reinforced  concrete  structures 
but  would  not  be  applicable  to  very  brittle  material  such 
as  unreinforced  concrete.  The  main  underlying  assumption 
is  that  the  structure  will  take  on  the  shape  as  proposed 
and  portions  of  the  structure  respond  almost  as  rigid  body 


motion  with  localized  failures  only  at  hinges.  Experi- 
1 9 

ments  ’ show  that  as  a small  explosive  device  is  moved 
closer  and  closer  to  a much  larger  plate  or  beam  localized 
failure  of  the  concrete,  i.e.,  actual  separation  of  concrete 
from  reinforcing  elements,  occurs  and  the* structure  simply 
does  not  contain  the  pressure.  This  type  failure  is  a 
mechanism  4 type  as  described  in  Section  2.1.  The  onset  of 
this  type  of  failure  appears  to  be  dependent  on  the  com- 
pressive strength  of  the  concrete  or  the  dynamic  compressive 
pressure  that  the  concrete  can  support  before  extreme  crack- 
ing occurs.  As  a general  rule  an  approximate  dynamic  load 
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factor  of  1.5  to  2.0  is  usually  applied  for  many  preliminary 
design  problems  and  this  type  of  thinking  appear:;  reasonable 
in  this  study.  Based  on  this  the  upper  bound  of  the  maximum 
pressure  of  the  spatial  variations  (P^+Pj.)  should  be  set  at 
approximately  1.5  times  the  compressive  strength  of  the  con- 
crete. This  rule  is  by  no  means  meant  to  be  a hard  and  fast 

rule  but  appears  to  be  reasonable  based  on  experimental 
• 1 9 

observations  ’ of  both  damage  and  pressure  time  histories. 

The  limit  of  rotation  or  center  point  defection  is 
based  on  the  failure  criterion  discussed  in  Section  2.3.4. 
Large  rotations  displacements  and  strains  occur  and  were 
accommodated  in  the  basic  assumptions.  In  the  strain  to 
fracture  calculations  the  expression  6z/2  was  considered 
to  be  small  in  comparison  to  unity.  This  more  than  likely 
should  limit  the  rotation  to  something  less  than  a third 
of  a radian. 

3.2  Numerical  Analysis  Procedure 

The  general  equations  as  developed  in  Section  2.3 
are  rather  complex  and  closed  form,  solutions  could  be  ob- 
tained given  time.  In  fact  closed  form  solutions  are  given 
for  uniform  loading  of  metal  1 -ame1'  and  parabolic  loadings 
of  concrete  b-ams  . However,  for  the  nonlinear  blast  load- 
ing of  plates  the  closed  form  solutions  become  very  cumber- 
some and  due  to  discontinuities  1;  the  solutions  cue  io 
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assumed  time  functions  and  geometric  response  shape  a numeri- 
cal technique  was  developed  to  solve  both  plate  and  beam 
problems.  The  steps  required  for  the  solution  are  as  follows: 

1.  Specify  beam  or  plate  problem. 

2.  Compilation  of  basic  data  on  idealized  structure: 
geometric  size,  pressure  loadings,  mechanical 
properties,  position  of  load.  (These  become 
input  parameters . ) 

3.  Determine  final  hinge  positions.  Tor  beams  +his 

is  automatically  the  half  span  of  the  beam,  and 
for  plates  this  means  solution  of  Equation  (26) 
or  (27)  for  z.  The  value  of  z determined  then 
becomes  a constant  for  the  problem. 

4.  Determine  the  initial  position  of  the  hinge 
location  as  an  indicator  of  which  equation  of 
motions  are  to  be  used.  The  initial  condition 
of  Equation  (14)  is  used  for  this  portion  of  the 
solution.  For  the  beams  Equation  (14)  is  applied 
to  Equations  (10)  and  (12)  or  Equations  (16).  For 
the  plates  the  initial  condition  is  applied  to 
Equations  (29)  and  (31)  or  Equations  (33)  and  (34). 

5.  From  No.  4 above  if  X’n  = a,(X=l)  then  response 
is  determined  from  a 6 equation  only.  For  beams 
the  response  is  based  on  Equation  (10)  or  (16a). 
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For  plates  the  response  is  obtained  using  Equa- 
tions (29)  or  (31).  If  0<x^<a,  (0<X<1)  then  the 
initial  response  is  governed  by  both  the  G and  6 
equations.  The  response  for  beams  is  then  deter- 
mined from  Equations  (10)  and  (12)  or  Equations 
(16).  For  plates  Equation  (29)  and  (31)  or  Equa- 
tions (33)  and  (34)  will  be  used  for  0<x^<a.  If 
the  initial  response  is  0<xp,<a  then  the  hinge  loca- 
tion moves  toward  the  center  position  and  when 
x^  = a the  solution  reverts  to  the  mechanism  1 mode 
of  using  only  the  G equations. 

6.  In  No.  5 above  the  proper  response  equations  of 
motion  are  selected  along  with  the  correspoinding 

W F,  WP,  PE  equations  and  they  are  solved  simultaneously. 

7.  Calculate  9 , 6 , 6 , WF, WP, PE, x^ , 6U . Print  data  for  cur- 
rent time,  compare  6 to  0U  and  print  asterisk  if  6>0U. 

8.  Return  to  next  case. 

The  numerical  technique  used  in  solving  for  the  initial 
hinge  location  x^  for  beams  and  plates  and  for  the  final  hinge 
location  zb  for  plates  is  a bisection  (binary  search)  method-^ 
described  in  some  detail  in  Section  A. 5.1  of  the  appendix. 

Since  the  rate  of  plastic  work,  rate  of  pressure  work  and 
rate  of  potential  energy  change  are  all  functions  of  0 and 
must  be  integrated  with  respect  to  time  numerically  they  are 

4G . 


combined  with  the  equations  of  motion  and  solved  simul- 
taneously using  a 4th  order  Rung  -Kutta  technique^1-.  This 
method  is  also  discussed  in  detail  in  Section  A. 5.2  of  the 
appendix . 


3.3  Numerical  Examples 
3.3.1  Introduction 

The  numerical  examples  included  here  are  response  cal- 
culations of  two  test  items  of  Reference  9.  The  test  items 
as  shown  typically  in  Figure  17  are  completely  buried  with 
soil  before  testing.  The  explosive  is  placed  so  as  to  be  on 
a line  normal  to  and  at  the  center  of  the  test  face.  The 
test  items  are  essentially  six  sided  boxes  of  36x36xl6in 
( . 91x . 91x . 41m)  and  48x36xl6in  (1. 22x. 91x. 41m)  with  the  side 
opposite  the  test  face  open.  The  side  walls  and  the  test 
face  are  4 in  (.lm)  thick.  For  numerical  analysis  the  test 
face  was  modelled  as  a 34x34x4in  ( . 9 lx . 9 lx. lm)  plate  for 
test  case  1 and  44x34x4in  ( 1 . 12x . 9 lx . lm)  for  test  case  2 
with  1%  tensile  reinforcement  ratio  for  both.  Reinforcing 
elements  are  spaced  2 in  (5.1cm)  apart  through  the  thickness 
and  l.Oin  (2.5cm)  from  the  plate  surface  giving  a reinforced 
distance  d (see  Figure  3)  of  3. Gin  (7.6cm)  for  each  slab. 

Based  on  pressure  measurements  for  tests  of  References 
1 and  9 a linear  pressure  (3=0)  of 
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Figure  17.  Schematic  of  test  specimen  used  in 
numerical  example. 


P(x)  = 2100  + 4300-  psi 

<3. 

P(y ) = 2100  + 4300^-  psi 

cL 


(4 


was  used  for  the  spatial  loading  of  each  slab.  A general 
time  function,  with  a = 0.  of 


was  used  for  the  pressure  time  variation.  The  general  shapes 
of  the  spatial  variation  and  time  variation  are  shown  in  Figure 
18.  From  Reference  9 the  concrete  compressive  strength  was 


time,  millsec 
b)  time  variation 


Figure  18.  Pressure  loading  for  numerical  example. 


6000  psi  (41.4MPa)  and  the  reinforcement  strength  v/as  70,000 
psi  (483.6  MPa).  The  time  step  (TINCR)  for  the  numerical 
analysis  was  chosen  as  0.5xl0~4  sec  with  a data  print 
frequency  (TPRINT)  every  0.2xl0-^  sec  and  a maximum  real 
time  response  (TMAX)  of  0.1  sec. 


3.3.2  Results 

The  results  of  the  two  numerical  examples  are  shown  in 
Tables  I and  II.  For  the  tabulated  results  of  TIME,  THETA, 
etc.,  the  time  is  given  only  in  the  upper  left  hand  column. 
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The  rows  in  the  lower  columns  PRESSURE  WORK,  PLASTIC  WORK, 


ere.,  correspond  to  the  rows  of  the  TIME  column.  The  pressure 
loadings  used  for  the  numerical  examples  were  taken  from 
References  1 and  9.  All  the  pressure  data  for  References  1 
and  9 were  collected  and  using  a regression  analysis  two 
regions  of  pressure  versus  standoff  distance  were  established. 
The  resulting  curves  are  shown  in  Figure  19. 

3.3.3  Discussion 

The  final  deformed  shapes  of  the  slabs  are  shown  in 
Figure  20.  It  appears  that  the  blast  load  in  each  case  was 
slightly  above  the  center  of  the  test  face  and  is  typical  of 
explosions  in  disturbed  soil  where  the  explosive  bubble  tends 
to  rise  into  the  top  loose  soil.  The  analysis  showed  the 
initial  hinge  location  vxas  approximately  11.5  in  (29.2cm) 
from  the  edge  for  each  slab.  In  a cross  section  across  the 
center  line  of  the  slab,  the  deformed  shape  after  some  dis- 
placement has  occurred  would  be  similar  to  the  sketch  of 
Figure  21.  The  rotation  occurs  at  A and  C of  Figure  21  and 
tensile  cracks  form  in  the  concrete  at  and  near  A and  com- 
pression cracks  and  dislodged  material  occurs  at  the  hinge 
at  C . In  general  this  appears  to  have  occurred  in  both  test 
specimens  of  Figure  20.  The  analysis  predicts  a maximum 
deflection  of  4 . 4in  (11.2cm)  with  a maximum  experimental 
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TABLE  I. 


RESULTS  UE  NUMERICAL  EXAMPLE  NO.  1 


TEST  CASE  i 


CALCULATIONS  uN  A CONCRETE  PLATE 
34  X 34  X 4 LINEAR 


INPUT  VALUES 


PLATE  HALF  WIDTH  CR  BE  A 1 HALF  SPAN,  IN. 
BEAH  CR  PLATE  THICKNESS,  IN. 

LENGTH  TO  WIDTH  RATIO,  DI ME  NS i 0 NL E SS 
MASS  PER  UNIT  AREA,  L 6S-S EC . SCO /I N . CU B ED 


(A)  17000000 

(H)  4.0000 0C0 

( A F ) 1.0GGQ00C 

(XM)  . 6996G0UG E-03 


MAXIMUM  DISTRIBUTED  PRESSURE  LOAD,  PS1 
UNIFORM  PRESSURE  LOAD,  PS  I. 

PRESSURE  DECAY,  DIMENSIONLESS 
PRESSURE  DURATION,  SEC. 


(PC) 
(FE) 
( ALPHA) 
(TAU) 


4300.0000 

2100.0000 

0 . 

. 5UC000GGE-Q 3 


SUPPORT  FACTOR!  1=SIMPLY,  2= CLAMPED  (F) 
WAVE  FUNCTION:  1=GENERAL,  2=SQJAKE  (WAVEFN) 
WEIGHT  VECTOR:  0=VERT,  1=EXP  BL  W,  -1=EXP  A3  V (AIDA) 
SPATIAL  PRESSURE  DECAY  CONSTANT .DIMENSION LESS (BETA) 


2.00JJO0C 
1. OGG0GOO 

0 . 

0 . 


CONCRETE  COMPRESSIVE  STRENGTH,  PS  I . (SIGMAC) 
REINFORCED  STEEL  YIElD  STRESS,  FSI.  (SIGMAh) 
REINFORCEMENT  RATIO  IN  TENSION,  DIMENSIONLESS  (O) 
REINFORCING  DISTANCE,  IN.  ( D) 


LOGO. 0000 
7 000J , 000 
.10000000 E-Dl 
3.0000000 


TIME  INCREMENT,  SEC.  (TINCR) 
TIME  MAXIMUM,  SEC.  ( T M AX ) 
TIME  STEF  INTERVAL  PER  PRINTED  LINE,  SEC.  (TPRINT) 


. 50000 COCE -04 
.10000200 
. 20002000  E-03 


COMPUTED  CONSTANT  VALUES 


PLATE  HALF  LENCTH,  IN'.  (L) 
RATIO  UF  FINAL  HINGE  lOC  TO  B,  DIMENSIONLESS  (2) 
HINGE  MOMENT,  IN. -LBS. /IN.  (XMU) 
WEIGHT  PER  UNIT  AREA,  L8S./IN.SQ.  (K) 
ORIGINAL  HINGE  LOCATION,  IN.  ( X H ) 


17. 002000 

I. 000  0 2Cu 
5279.7150 
• 34760  544 

II. 642934 


CLAMPED-SUPPORTED  GENERAL  TIME  FUNCTION 


4 5 . 


VERTICAL  WAL 


LINEAR  LOAu 


TABLE  I.  RESULTS  NUMERICAL  EXAMPLE  NO.  1 (Concluded) 


CALCULAT1Q 

NS  ON  A CONCRETE 

PLATE 

TEST  CAS  F 

1 <4  X 34  X 

4 L 1 N E A R 

C L A M f 

- L- SUPPOF  TEu 

GlN  ERA  L TIME  FUNCTION 

TIME 

THETA 

MIuPT.  \l  EL. 

M1DFT.  OELT 

1 

(SECONDS  ) 

(RADIANS) 

(IN. /SEC.  ) 

(INCHES) 

0 • 

0. 

0. 

0. 

• 2 00  00  300E-03 

. 9 01'J  1178F-02 

976. 8 0 1 6 fa 

. 10  572560 

.LOO  0 0 0 00E -0  3 

. 2 4 95 7 7 6 HE- 0 1 

14o8.2165 

. 35  8 36  3 34 

.60000000  E -0  3 

.539ul49GE-ui 

153C. uUol 

• bb  2 30  3 37 

. 830000  0 0E-03 

. 76376661E-01 

1530. 0 0 61 

. 9G830469 

.luOGOOOTE-02 

. *1728277 4 1 — 0 1 

1530. OC'ol 

1.2743058 

. 12000000E-02 

• 11670349 

153  U . G 061 

1.6803070 

.140 00000 E-02 

• 13471441 

1530. Gufal 

1 .8863  062 

. 160  OOOGJE-G2 

. 1 61 38235 

153  0. G Gbl 

2.  1 9230  95 

. 180  00  COOE-Q  2 

. 166766  20 

1630.0061 

2.4983107 

. 200  00  OOOE-O  2 

. 1 80  51173 

163  0.  0 061 

2.  e043H9 

. 220  00UG3E- 02 

.19386114 

153U . C 061 

3.1103121 

.24000000E-02 

• 20  564225 

1 530  . 0 0 61 

3.4163l43 

. 260  00003E-02 

.21627202 

864.78575 

3.67  66  2 1 3 

.28000CG0E-02 

. 22 575492 

757.3  j 824 

3. 83  783  3 7 

. 30000000E -02 

. 2 3409103 

659.83073 

3.97954/6 

•32000000E-02 

. 24 1 280  35 

562.35322 

4. 1 01  7669 

. 340  OOOOOE-D2 

. 24732287 

464 .87571 

4.2044888 

. 360  C0G00E-CI2 

. 26221 650 

3t  7 . 3 9 8 2 C 

4.287/162 

.38JQO0QOE-02 

• 2 5 5 9 6 7 54 

269. 92069 

4. 3514481 

. 400  00  000E-02 

. 25  85  6 958 

172 .4431 7 

4.3956845 

. 420  OOOOOE-O  2 

. 26002602 

74. 965665 

4.4204254 

I PRESSURE  WCRK  PLASTIC  6QRK  KINETIC  ENERGY  HINGE  LOCATION 


(IN.  -LBS.  ) 

(IN.-LoS.) 

(IN.  — L u s • ) 

< INCHES) 

c . 

0. 

0 . 

11. 542934 

177758.94 

12952. 197 

16480b. 74 

11.  7223  88 

395472 .30 

43021.825 

352430 .47 

11. 962284 

427914.43 

77  40  6. 985 

3 50  50  7.  44 

12.  28/2  52 

427914.43 

1 0 96 83 . 18 

318231. 24 

12.676017 

427914.43 

I397u6.09 

28820o.34 

13.058987 

427914  .43 

167  595 .84 

260  51 8. 59 

13.541215 

427914.43 

193461.00 

2 3445  3.  43 

14.002276 

427914.43 

217397.54 

210516.89 

14. 48195u 

427914.43 

239488.59 

186  425.  84 

14. 981007 

427914.43 

259804. 17 

168110.26 

15.506955 

427914.43 

278400 .59 

14  3 613.  84 

16.044 U 25 

427914.43 

295319.24 

132596. 19 

lo. 6l29ui 

427914.43 

310  584.4? 

117229. 98 

17. QOUOUd 

42791  4.43 

324202.69 

103711. 74 

l/.OuOuGU 

427914.43 

336174.03 

91740.397 

l/.OOtOjl 

42791 4. 43 

346498.43 

81415.945 

17. QGG0G0 

427914.44 

355176. 04 

72738. 383 

17.000000 

427914.43 

3622C6. 72 

657  07  . / 13 

17.000000 

42791  4 .43 

367590.49 

60323.933 

17.0000CO 

427914.43 

371327.38 

565  87 . 0 45 

17.000000 

427914.43 

373417.38 

6 4 4 9 7 .048 

17. 0 0 0 0 0 0 

MAXIMUM  LE- 

LECTION = 4.4256708  AT  TIME  = 

.4-»GGGG0CE-Q2 

AN  ASTc-.lSK 

INDICATES  THAT  A 

REINFORCING  ELEMENT 

HAS  FRACTURED 

i 
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TABLE  II. 


RESULTS  OF  NUMERICAL  EXAMPLE  NO.  ? 


CALCULATIONS  CN  A C UN  CRETE  PLATE 
TEST  CASE  2 44  A 34  X 4 LINEAR 


INPUT  VALUES 


PLATE  HALF  WIDTH  CR  BEAM  HALF  SPAN,  IN. 
BEAM  OR  PLATE  THICKNESS,  IN. 

LENGTH  TO  WIDTH  RATIO,  DIMENSIONLESS 
MASS  PER  UNI  1 AREA,  L oS -S EC . S GJ / 1 N . C U6ED 


(A)  1 7 . G 0 0 0 G G 

(H)  4.GGJOOGC 
(AK)  1.2  90  0 GOO 
(XM)  . B99oG DGOE-O 3 


MAXIMUM  DISTRIBUTED  PRESSURE  LOAJ,  PS1. 
UN  IF  Jr-  M PRESSURE  LuAL,  P5I. 

PRESSURE  DECAY.  DIMENSIONLESS 
PRESSURE  DURATION,  SEC. 


( PC) 
(PE  ) 
( ALPHA) 
(TAU) 


43GC.UGDG 
21G  0. GGGO 

0 . 

. 5UC0CQ00E-Q  3 


SUPPORT  FACTOR*  1=SIMPLY,  2=CLAMPE0  <F)  2.GGJ0G00 
WAVE  FUNCTION:  1=GENERAL,  2=SQUAKE  (WAVEFN)  1.GJG00G0 
WEIGHT  VECTOR:  0=VEkT,  1=EXP  GlW,  -1=EXP  ABV  (AIDA)  u . 

SPATIAL  PRESSURE  DECAY  CONST ANT ,D 1 MENSICNLESS ( BET  A)  o. 


CONCRETE  COMPRESSIVE  SlKlNGTH,  F S 1 . (SIGMAC) 
REINFORCED  STEEl  YIELD  STRESS,  PS  I . (S1GMAR) 
REINFORCEMENT  RATIO  IN  TENSION,  DIMENSIONLESS  ( Q ) 
REINFORCING  DISTANCE,  IN.  (C) 


BOO  G. OGGG 
7GGG...GGG 
. 10 GJOG ODE-01 
3. OCUGOCG 


TIME  INCREMENT,  SEC.  ( T I NCR ) 
TIME  MAXIMUM,  SEC.  ( T M A X ) 
TIME  STEP  INTERVAL  PER  PRINTED  LINE,  SEC.  (TPRINT ) 


. 5GCjOOOOE-G4 
.10GG0GG0 
. 2UG JGUOQ  E-Q  3 


COMPUTED  CONSTANT  VALUES 


PLATE  HALF  LENGTH,  IN.  (G) 
RATIO  OF  FINAL  HINGE  LOC  TO  3,  DIMENSIONLESS  (Z) 
HINGE  MOMENT,  IN. -LOS. /IN.  (XMU) 
WEIGHT  PER  UNIT  AREA,  LBS./IN.SQ,  (W> 
ORIGINAL  HINGE  LOCATION,  IN.  (XH) 


21.930000 
. 87892855 
5279.7150 
• 347 60  544 
11. 4ob 1 99 


CLAMPED-SUFFORTEO  GENERAL  TIME  FUNCTION 


VERTICAL  WALL 


LINEAR  LOAD 


TABLE  II.  RESULTS  OF  NUMERICAL  EXAMPLE  NO.  2 (Concluded) 

CALCULATIONS  ON  A CONCRETE  PLATE 
44  X 34  X 4 LINEAR 


TEST  CASE  2 

clam ped-sup pot t eo 


general  time  function 


r i me 

(SECONDS) 


0 . 


• 2GG  00  00 0 E -0  3 
. 40000000--03 
.600  0000JE-03 
. 8D0  00O0  TE-03 
. lflt  00  00  9 E-0  2 
. 120  JO  0 0 3 E -0  2 
. 1400000 JE -02 
. 16000  OQ9E-Q2 
.18000000E-02 
. 200  OGGO0E-O2 
. 220000 J3E-02 
. 24000000E-02 

• 260  OJOOOE-02 
.280  Q0G00E-G2 
.3000000JE-02 
.3  20  0000 0E-02 
. 34000003E -02 
. 3600000 JE-02 
. 380  00  000E-02 
. 40  000000E-C2 
.420  03  0 JOE -02 
. 440 00 000 E -02 
. 460  00  uOOt -0  2 
•4800000QE-G2 
. 5 0 0 0 0 0 0 0 E - 0 2 
. 520  0000 JE-0  2 


THETA 

(RADIANS) 


. 88594534t-u  2 
. 29497509E-01 
. 53255951E-U 1 
. 757H1298E-ol 
. 96  96  25  9 3t  -b  1 
. 11686975 
. 1 3556763 
. 15311561 
. lo956757 
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AN  ASTERISK  INDICATES  THAT  A REINFORCING  ELEMENT  HAS  FRACTUktU 
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Figure  21.  Schematic  of  plate  cross  section  through 
middle  of  plate. 

deflection  of  2.8in  (7.11cm)  for  test  case  1.  For  test 
case  2 the  analysis  shows  a maximum  displacement  of  5.25in 
(13.3cm)  and  the  experimentally  observed  maximum  deflection 
was  approximately  8in  (20.3cm)  at  the  center.  The  analysis 
indicates  failure  of  a reinforcing  element  occurred  in  test 
case  2 and  not  of  1 and  Figure  20  shows  evidence  of  severe 
damage  of  test  case  2 and  not  of  test  case  1.  Probably  the 
most  important  result  of  these  analyses  is  that  the  analysis 
distinguishes  between  the  severity  of  damage  with  just  changes 
in  the  aspect  ratio  AR  of  the  slab. 


I 

I 

The  failure  criterion  as  giver  in  the  study  uses  an 
N ultimate  strain  eu  of  &2 based  on  a reported  elongation  range 

of  0.15  to  0.25  for  various  steels  . The  value  of  eu  was 
fixed  in  the  computer  program  as  0.2  and  appears  to  be 
justified  on  the  fact  that  the  major  portion  of  reinforced 
concrete  has  steel  as  reinforcing  elements. 

IV.  CONCLUSIONS  AND  RECOMMENDATIONS 

Experiments  using  small  explosive  devices  verify  that 
reinforced  concrete  plates  and  beams  exhibit  a plastic  hinge 
response  when  subjected  to  dynamic  loadings  of  such  devices. 

Using  failure  mechanisms  based  on  plastic  hinge  response, 
equations  of  motion,  for  reinforced  concrete  beams  and  plates, 
are  derivable  and  their  solutions  predict  responses  that  are 
in  reasonable  agreement  with  experiments.  These  solutions 
are  limited  to  simply  supported  and  fixed  edge  conditions  and 
their  accuracy  is  heavily  dependent  on  the  accuracy  of  the 
pressure  time  loading  function.  The  basic  assumptions  of 
these  analyses,  i.e.,  plastic  hinge  response,  neglecting 
elastic  motion  and  tensile  concrete  strength,  coupled  with 
the  variations  found  in  pressure  loading  predictions  could 
produce  overall  response  errors  of  j^50i. 

In  addition  application  of  the  derived  equations  of 
motion  are  limited  to  pressure  loadings  between  the  pressure 
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required  for  static  collapse  and  pressures  where  the  concrete 
begins  to  break  up  locally  before  the  plastic  hinge  response 
begins . 

It  is  recommended  that  the  study  of  blast  response  of 
reinforced  concrete  structures  be  continued  to  include  non- 
symmetric  loadings,  response  to  very  localized  or  very  near 
field  explosives  and  complete  shearing  of  small  structural 
elements  at  boundaries  (Mechanisms  3 and  4 of  Section  2.1). 
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APPENDIX 


COMPUTER  PROGRAM  SPECIFICATION  FOR  CONCRE 

A.l  Introduction 

The  following  describes  the  computer  program  named 
CONCRE,  Reinforced  Concrete  Beam  or  Plate  Center  Deflection. 
This  program  was  developed  and  used  to  test  the  theory  ex- 
plained in  the  earlier  part  of  this  document.  For  a listing, 
see  Section  A. 6. 

The  program,  CONCRE,  computes  the  deflection  of  the 
center  point  and  edge  rotation  of  a reinforced  concrete 
beam  or  plate  caused  by  a linear  or  blast  load.  The  main 
functions  of  the  program  are  DRIVER,  NEXT-CASE,  TIME-CONTROL, 
TIME-ZERO,  BISECTION,  TIME- STEP,  RUNGE-KUTTA,  CHECK-HINGE,  and 
PRINT-ROUTINE.  See  Diagram  A. 2.1  for  a description  of  the 
interaction  of  these  functions. 

This  program  was  developed  using  the  FORTRAN  Extended 
Compiler  and  the  N0S/BE1  Operating  System  on  the  CDC  6600 
Computer  at  Eglin  Air  Force  Base,  Florida.  It  was  written 
in  a top-down,  structured  manner,  containing  numerous  comment 
It  uses  Standard  FORTRAN,  with  the  exception  of  the  use  of 
asterisks  (*)  to  define  Hollerith  data  within  FORMAT  state- 
ments and  calls  to  certain  system  library  functions.  The 
program  is  complete  in  itself. 
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The  input  to  the  program  should  be  in  eighty-column 
card  format.  The  only  output  is  a printed  listing,  which 
includes  the  input  as  read,  the  results  as  requested  and 
computed,  and  error  messages,  if  any.  For  a further  descrip- 
tion of  the  input,  see  Section  A. 4.1.  Section  A. 4. 2 further 
describes  the  output.  A sample  run  complete  with  input  and 
output  is  given  in  Section  3.3. 

A. 2 Functional  Description 

The  nine  main  functions  of  the  program,  CONCRE,  are 
described  in  the  subparagraphs  below.  The  interaction  be- 
tween these  functions  is  shown  in  Figure  A. 2.1.  An  abbre- 
viated name  is  given  for  each  function  for  clarity  and  simpli- 
fication in  the  diagrams.  In  most  cases,  these  abbreviated 
names  match  the  name  of  one  of  the  program's  subprocedures, 
although  a function  may  be  made  up  of  more  than  one  procedure. 
Figure  A.  2. 2 show’s  the  relationship  between  the  21  subprocedures. 
For  further  information  on  each  function,  see  the  comments 
embedded  in  the  listing  (Section  A. 6)  for  each  procedure 
within  the  function. 

A . 2 . 1 Driver 

The  function,  DRIVER,  drives  the  program  taking  care  of 
initialization  and  clean  up  activities  as  well  as  controlling 
the  main  flow  of  the  program.  This  is  the  main  function  of 
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the  program  and  begins  when  the  program  begins  execution. 

It  is  terminated  when  an  End-of-File  condition  is  returned 
by  function  NEXT-CASE. 

This  function  is  made  up  of  the  two  FORTRAN  procedures, 
CONCRE  and  DRIVER.  It  calls  the  functions,  NEXT-CASE  and 
TIME-CONTROL.  Flow  Diagrams  of  the  two  procedures,  CONCRE 
and  DRIVER,  are  shown  in  Figures  A. 2. 1.1  and  A. 2.1. 2. 

A. 2. 2 Next-Case 

The  function,  NEXT-CASE  (or  NXTCAS) , finds  the  next  case 
in  the  stream  of  input,  reads,  stores,  and  prints  out  the 
data  associated  with  that  case.  If  errors  are  detected  in 
the  input  data,  proper  error  messages  are  printed,  and  the 
Input-Error-Flag  is  returned  to  the  DRIVER  function  as  set 
so  that  this  case  will  be  terminated  and  NEXT -CASE  recalled. 
Should  the  End-of-File  be  found  before  finding  the  first 
card  of  the  next  case,  the  function  returr s control  to  the 
DRIVER  function  for  a normal  termination.  Should  the  End- 
of-File  be  found  before  reading  the  last  data  card  of  the 
next  case,  error  messages  will  be  printed  before  return  to 
the  DRIVER  function. 

This  function,  NEXT-CASE,  is  called  by  and  returns  to  the 
DRIVER  function.  It  sets  the  Input-Error-Flag  and  End-of- 
File-Flag  used  by  the  DRIVER  function.  It  calls  no  other 
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It  is  made  up  of  the  following  four  FORTRAN 
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procedures : 

NXTCAS , which  drives  the  function; 

GETITL,  which  searches  for  the  first  card  of  the  case; 

RWDATA,  which  reads  and  prints  the  input  data;  and 

PAGE,  which  prints  the  headings  on  the  output. 

Flow  diagrams  for  the  first  three  procedures  are  shown  in 
Figures  A. 2. 2.1  through  A. 2. 2.3.  Figure  A. 2.9.2  provides 
the  Flow  Diagram  for  the  PAGE  procedure. 

A. 2. 3 Time-Control 

The  function,  TIME-CONTROL  (or  TCNTPL),  controls  the 
production  of  computed  results  for  each  case  starting  with 
initialization  at  time  zero  through  all  the  integration 
time  steps  until  time-max  is  reached  or  the  maximum  deflec- 
tion is  found  (i.e.,  the  velocity  becomes  negative),  which- 
ever happens  first.  If  problems  arise,  such  as  a bad  hinge 
location  or  negative  constants  being  cc  ted,  appropriate 
error  messages  are  written  and  co.  ; ol  is  returned  to  the 
calling  function. 

This  function,  TIME-CONTROL,  is  called  by  and  returns 
to  the  DRIVER  function.  It  sets  the  Case- is-Done-Flag  if 
computation  is  completed,  if  disasterous  errors  are  found, 
or  if  the  case  has  been  tried  a maximum  number  of  times. 
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Figure  A.  2. 2.1.  Flow  Diagram  of  Subroutine  NXTCAS 


GE7ITL 


Figure  A. 2. 2. 2.  Flow  Diagram  of  Subroutine  GETITL 


Figure  A. 2. 2. 3.  Flow  Diagram  of  Subroutine  RWDATA 


The  parameter,  NTRIES,  which  is  returned  to  the  calling 
function,  counts  the  number  of  times  the  current  case  has 
called  this  function,  TIME-CONTROL. 

Two  other  functions  are  called  by  the  TIME-CONTROL 
function;  TIME-ZERO,  which  handles  initialization  and  re- 
sults at  time  zero;  and  TIME-STEP,  which  handles  the  computa- 
tions for  the  rest  of  the  integration  time  steps  , t] , . . . , 

tmax*  The  parameter,  NTRIES,  is  passed  to  (but  not  altered 
by)  the  function,  TIME-ZERO,  to  eliminate  re-initializing 
case  constants  for  the  second  or  succeeding  tries.  The 
Loop-Control-Elag  is  returned  as  set  by  the  TIME-STEP 
function  if  the  loop  should  be  terminated  for  any  reason. 

Only  one  FORTRAN  procedure  makes  up  this  function; 
namely,  TCNTRL.  See  Figure  A. 2. 3.1  for  a Flow  Diagram  of 
this  procedure. 


A. 2.4  Time-Zero 

The  function,  TIME-ZERO  (or  TZERO),  initializes  all 

variables  and  constants  to  be  used  by  the  case,  including 

finding  the  value  of  the  original  hinge  location,  xh  , and 

0 

for  a plate  case,  z,  the  ratio  of  the  final  hinge  location 
to  plate  half  length.  It  also  provides  that  the  headings 
and  first  line  of  output  results  for  the  case  be  printed. 


Should  invalid  constants  (i.e.,  negative)  be  calculated,  or. 


rCNTRL(NTRIES) 


Figure  A. 2.3.1.  Flow  Diagram  of  Subroutine  TCNTRL 


for  a plate  case,  a value  of  2 not  converge,  the  Case-is- 


Done-Flag  will  be  set. 


Called  by  the  function,  TIME-CONTROL,  this  function, 


uses  the  parameter,  NTRIES,  to  prevent  recalculating  constants 


or  the  original  hinge  location  after  the  first  try  for  a 


given  case.  This  parameter  is  not  altered.  The  setting  of 


the  Case-is-Done-Flag  as  above  assures  the  case  will  be  ter- 


minated upon  return  to  the  calling  function,  if  necessary. 


Three  other  functions  are  called  by  the  TIME-ZERO  func- 


tion: BISECTION,  a root-finding  function  used  to  find  the 


values  of  x-  and  z;  CHECK-HINGE,  to  assure  the  proper 


Mechanism-Flag  setting  for  the  newly-found  hinge  location, 


and  PRINT-ROUTINE,  to  print  the  headings  and  first  line  of 


output  for  the  given  case.  The  Non-Convergent-Flag  will  be 


returned  by  the  BISECTION  function  if  Xn  or  z would  not 


converge.  The  name  of  the  proper  procedure  to  be  used  by 


the  BISECTION  function  in  determining  Xn  or  z is  passed  to 


that  function  as  an  argument  (which  is  not  altered). 


Five  FORTRAN  procedures  make  up  the  TIME-ZERO  function: 


TZERO,  for  general  initialization  and  to  drive  the 


other  procedures; 


BTZERO,  for  initializing  beam  constants  and  variables; 


PTZERO,  for  initializing  plate  constants  and  variables; 


CALXKO,  for  deriving  the  original  hinge  location;  and 


CALBAR,  to  assist  in  initializing  .!»<•  constants 
for  mechanism  1 Runge-Kutta  computations. 

Flow  Diagrams  of  these  first  four  procedures  will  be  found 
in  Figures  A.  2.4.1  through.  A.  2. 4. 4.  The  Flow  Diagram  for 
the  procedure,  CALBAR,  will  be  found  in  Figure  A. 2. 7. 3. 

A. 2. 5 Bisection 

The  function,  BISECTION  (or  BISECT)  is  a root-finding 
technique.  Using  the  BISECTION  (or  Binary  Search)  Method, 
this  function  derives  the  root  of  one  of  three  different 
formulas  to  find  the  values  of  Xq  and  z.  If  the  root  cannot 
be  found  or  does  not  converge,  the  Error-Flag  (Non-Convergent- 
Flag)  is  returned  to  the  calling  function. 

Called  by  the  TIME-ZERO  function,  this  function  uses 
parameters  to  define  the  likely  interval  within  which  the 
root  should  be  and  the  formula  to  be  used.  It  returns  the 
root,  if  found,  and  the  Error-Flag,  if  set,  to  the  calling 
function.  The  BISECTION  function  calls  no  other  function. 

This  function  consists  of  five  FORTRAN  procedures: 

BISECT,  the  actual  root-finding  procedure,  which 
can  stand  alone; 

FTNZ , the  procedure  defining  the  z formula  for  a 
Plate  Case; 

BFTNX,  the  procedure  defining  the  Xq  formula  for  a 
Beam  Case; 
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Figure  A. 2. 4. 2.  Flow  Diagram  of  Subroutine  BTZEKO 
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Figure  A. 2. 4. 4.  Flow  Diagram  of  Subroutine  CALXHO 
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PFTNX,  the  procedure  defining  the  Xq  formula  for  a 
Plate  Case;  and 

CALBAR,  used  by  PFTNX  to  calculate  0 and  6 at  time 
zero . 

The  Flow  Diagrams  for  the  first  four  procedures  are  found 
in  Figures  A. 2. 5.1  through  A. 2. 5. 4.  The  Flow  Diagram  for 
CALBAR  is  Figure  A. 2.7. 3. 

A. 2 . 6 Time-Step 

The  function,  TIME-STEP  (or  TSTEP),  controls  the  compu- 
tation of  results  for  a single  integration  time  step.  It 
drives  the  other  functions  needed  to  compute,  check,  and 
print  these  results. 

This  function,  TIME-STEP,  is  called  by  the  TIME-CONTROL 
function.  The  Loop-Control-Flag,  LOOPFG,  is  a parameter  to 
this  function  and  is  returned  as  set  should  the  loop  be 
terminated  for  any  reason.  The  Case-is-Done-Flag  will  also 
be  set  when  the  last  time  step  needed  is  completed. 

Three  other  functions  are  called  by  the  TIME-STEP 
function:  RUNGE-KUTTA,  to  compute  the  equations  of  motion 

and  work  for  the  current  time  step;  CHECK-HINGE,  to  compute 
and  check  the  Mechanism-Flag  setting  of  the  new  resultant 
hinge  location;  and  PRINT-ROUTINE,  to  print  out  the  results 
for  this  time  step,  if  requested.  The  function  CHECK-HINGE 


BISECT ( GUESS1 , GUESS2 , ROOT , ERRFLG , FTN ) 


Figure  A. 2. 5.1.  Flow  Diagram  at  Subroutine  BISECT 


will  only  be  used  while  the  case  remains  in  Mechanism  2. 

All  the  results  are  passed  in  COMMON  blocks. 

The  TIME-STEP  function  is  made  up  of  only  one  FORTRAN 
procedure,  TSTEP.  Figure  A. 2. 6.1  shows  the  Flow  Diagram 
of  this  procedure. 

A. 2. 7 Runge-Kutta 

The  function,  RUNGE-KUTTA  (or  RUNGEK) , uses  the  fourth- 
order  Runge-Kutta  technique  to  solve  the  five  simultaneous 
differential  equations  of  motion  and  work  (two  of  which  are 
second-order  differential  equations)  for  the  given  time-step. 

A more  detailed  explanation  of  the  use  of  this  technique 
appears  in  Section  A. 5.2. 

This  function  is  called  by  the  TIME-STEP  function  and 
returns  computed  results  through  the  COMMON  block,  RESULT. 

It  calls  no  other  function. 

Three  FORTRAN  procedures  make  up  the  RUNGE-KUTTA  function 
RUNGEK,  which  drives  the  computation  procedure  four 
times  and  puts  together  the  resultant  pieces 
according  to  the  Runge-Kutta  formula; 

COMP,  which  actually  performs  the  computations  of 
the  five  differential  equations  using  parameters 
provided  by  RUNGEK;  and 

CALBAR,  which  provides  pieces  of  the  0 and  6 calcula- 


tion for  plate  cases. 


TSTEP(LOOPFG) 


Figure  A. 2. 6.1.  Flow  Diagram  of  Subroutine  TSTEP 


The  Flow  Diagrams  for  these  three  procedures  are  shown  in 


Figures  A. 2.7.1  through  A. 2. 7. 3. 


■ 

y 
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A. 2. 8 Check-Hinge 

The  function,  CHECK-HINGE  (or  CHEKXH)  computes  the  next 
value  of  the  hinge  location  for  a case  in  Mechanism  2 or  at 
time  zero.  If  the  hinge  has  moved  within  2%  of  the  center 
of  the  plate  or  beam,  the  case  is  considered  to  be  in  Mech- 
anism 1. 

This  function  is  called  by  both  the  TIME-ZERO  and  TIME- 
STEP  functions.  It  changes  the  setting  of  the  Mechanism- 
Flag  from  2 to  1,  if  the  hinge  is  within  2%  of  the  center. 

If  the  hinge  is  negative,  the  Case-is-Done-Flag  is  set  as 
the  data  must  be  bad.  The  Bad-Hinge-Flag  is  set  if  the 

S' 

hinge  is  negative  or  goes  beyond  2%  of  the  center  of  the 
beam  or  plate.  In  the  latter  case,  the  time  step  increment 
is  halved  so  that  the  next  try  for  this  case  will  have  a 
smaller  time  step,  hopefully  insuring  better  results. 

The  CHECK-HINGE  function  calls  no  other  functions  and 
consists  of  the  single  FORTRAN  procedure,  CHEKXH.  The  Flow 
Diagram  for  the  procedure  will  be  found  in  Figure  A. 2. 8.1. 

A. 2. 9 Print-Routine 

The  last  function,  PRINT-ROUTINE  (or  PRINTR) , provides 
for  the  printing  of  the  results  calculated.  It  is  called  by 
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Figure  A. 2. 7.1.  Flow  Diagram  of  Subroutine  RUNGEK 


Figure  A. 2. 7. 2.  Flow  Diagram  of  Subroutine  COMP 
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Figure  A. 2. 7. 3.  Flow  Diagram  of  Subroutine  CALBAR 
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both  the  TIME-ZERO  and  TIME-STEP  functions,  using  informa- 
tion passed  through  the  COMMON  blocks,  PRINTS  and  RESULT. 

It  calls  no  other  function. 

Two  FORTRAN  procedures  make  up  the  PRINT-ROUTINE 
function: . 

PRINTR,  which  prints  a single  line  of  results, 
adding  to  the  number  of  lines  printed;  and 
PAGE,  which  provides  headings  on  a new  page  of  out- 
put and  resets  the  line  count  whenever  a new 
page  is  needed. 

The  Flow  Diagram  for  these  procedures  can  be  found  in  Figures 
A. 2 . 9 . 1 and  A. 2 . 9 . 2 . 

A. 3 Storage  Allocation  , 

This  program  was  written  in  Control  Data  Corporation 
(CDC)  FORTRAN  Extended,  Version  4,  an  extension  of  ANSI 
FORTRAN  IV,  and  was  tested  on  the  N0S/BE1  operating  System 
of  the  CDC  6600  at  Eglin  Air  Force  Base  in  Florida.  Thus, 
all  figures  which  follow  refer  to  the  manner  in  which  storage 
was  allocated  on  that  machine  and  can  only  serve  as  a model 
for  those  wishing  to  use  the  program  elsewhere. 

All  data  base  tables  and  constants  are  internal  to  the 
program  and  are  described  in  detail  below.  A more  precise 
specification  of  storage  allocation  can  be  obtained  by  request- 
ing a Load  Map  during  the  loading  of  the  program. 
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Flow  Diagram  of  Subroutine  PAGE 


A. 3.1  Program  Variables 

All  the  variables  used  by  the  program  are  defined  in 
Table  A. 3. 1.1.  This  acts  as  a glossary.  The  variables 
are  listed  alphabetically  with  a brief  description  devoted 
to  each  one.  The  first  part  of  the  description  indicates 
the  sub-procedures  (in  parenthesis)  or  COMMON  block  (with- 
in slashes)  in  which  the  variable  is  used.  An  asterisk 
before  the  variable  name  indicates  that  it  is  an  input  var- 
iable. If  the  variable  is  actually  a table,  the  dimension 
is  provided  in  parenthesis  following  the  variable  name.  A 
separate  list  of  all  dimensioned  variables  appears  in 
Table  A. 3 . 1 . 2 . 

A. 3. 2 Common  Blocks 

Tables,  variables,  and  constants  that  are  used  by 
several  subprocedures  are  stored  in  labeled  COMMON  blocks. 
Table  A. 3. 2.1  outlines  the  data  flow  between  these  COMMON 
blocks  and  the  individual  subprocedures.  It  also  indicates 
the  length  of  each  COMMON  block. 

A brief  description  of  the  purpose  of  each  COMMON 
block  follows: 

FLAGS,  contains  the  flags  and  counts  used  to  con- 
trol the  flow  of  the  program; 

BARS,  contains  the  barred  elements  (A,...,U)  used  in 
computing  the  equations  of  motion  of  a plate  case; 
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TABLE  A. 3. 1.1.  LIST  OF  PROGRAM  VARIABLES 

*A  /INPUTS/  Beam  Half  Span  or  Plate  Half  Width, 

a,  in. 

AA(5)  (COMP)  Dummy  array  parameter  returning  first 

and/or  second  order  differential  equation 
values  to  Runge-Kutta  routine 

ACUBE  /CONSTS/  a3 

*AIDA  /INPUTS/  Weight  vector  for  plate  case,  n: 


0,  vertical  wall 

1,  horizontal  slab  with  explosive  below 

-1,  horizontal  slab  with  explosive  above 

i 

*ALPHA  /INPUTS/  Pressure  Decay  Constant,  a,  dimen- 


sionless 

*AR  /INPUTS/  Aspect  ratio,  length  to  width,  AR, 

dimensionless 

ARG ( 5 ) (RUNGEK  and  COMP)  Arguments  passed  by  Runge- 

Kutta  routine  to  computation  routine  in  order 
to  obtain  new  values  for  the  first  and  second 
order  differential  equations: 

ARG ( 1 ) = t 

J ARG (2)  = 6 

ARGO)  = 6 
ARGO)  = 6 
ARG ( 5 ) = 8 


1 


i 

l- 


TABLE  A. 3. 1.1.  LIST  OF  PROGRAM  VARIABLES,  Continued 


ARSQ 

ARSZP1 

ARZSP1 

ASQ 

ATHEDB 


ATHEDL 


ATHED1 


ATHED2 


AVDOT 


AWFDB 


/CONSTS/  AE2 
/CONSTS/  AR2Z  + 1 
/CONSTS/  AR  Z2  + 1 
/CONSTS/  a2 

/COMPS/  Constant  portion  of  first  part  of  0 
computation  for  a Beam  Case  with  Blast  load 
in  Mechanism  1 

/COMPS/  Constant  portion  of  first  part  of  8 
computation  for  a Beam  Case  with  Linear  load 
in  Mechanism  1 

/COMPS/  Constant  portion  of  first  part  of  0 
computation  for  Plate  case  in  Mechanism  1, 
Linear  or  Blast  load 

/COMPS/  Second  constant  portion  of  0 computa- 
tion for  Beam  or  Plate  case  in  Mechanism  1, 
Linear  or  Blast  load 

• • 

/COMPS/  Constant  portion  of  V(PE)  computation 
for  Beam  or  Plate  case  in  Mechanism  1 
/COMPS/  Constant  portion  of  Wf  computation 
for  Beam  or  Plate  Case  in  Mechanism  1,  Blast 


N 


yo . 


Lr. 


load 


TABLE  A. 3.1.1.  LIST  OF  PROGRAM  VARIABLES,  Continued 


i 

i 

i 

, 

I 


i 


AWFDL 

AO  (5) 

Al(  5) 

A2  ( 5 ) 
A3  ( 5 ) 

B 

BADXFG 

BARA 

BARA1 

BARB 

BARC 


/COMPS/  Constant  portion  of  WF  computation 
for  Beam  or  Plate  case  in  Mechanism  1,  Linear 
load 

(RUNGEK)  Computation  results  for  first  of  the 
four  Runge-Kutta  formulas,  aQ , for  all  five 
differential  equations 

(RUNGEK)  Computation  results  for  second  of  the 
four  Runge-Kutta  formulas,  a1 , for  all  five 
differential  equations 

(RUNGEK)  Third  Runge-Kutta  formula  results,  a2 
(RUNGEK)  Fourth  Runge-Kutta  formula  results,  ag 
/CONSTS/  Plate  Half  Length  or  Beam  Width,  b,  in. 

A 

/FLAGS/  Bad-Hinge-Flag  set  to  1 by  CHEKXH  if  a 
bad  hinge  location  is  computed;  otherwise,  zero 
(Integer) 

/BARS/  Used  to  compute  0 for  plate  case,  blast 
load , A 

/CONSTS/  Constant  portion  of  A for  given  plate 
case 

/BARS/  Used  to  compute  0 for  plate  case,  blast 
load , B 

/BARS/  Used  to  compute  0 for  plate  case,  C 
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TABLE  A. 3. 1.1.  LIST  OF  PROGRAM  VARIABLES,  Continued 


BARC1 

/CONSTS/ 

Constant  portion  of 

C for 

given 

plate 

case 

BARD 

/BARS/ 

Used  to  compute  0 for 

plate 

case , 

D 

BARDEL 

/BARS/ 

Used  to  compute  6 for 

plate 

case : 

1*  K, 

if  blast  load 

T*U , 

if  linear  load 

BARDNM 

/BARS/ 

Common  denominator  used  in  calculating 

A,  D,  H, 

and  R for  a plate  case  0 computation 

BARG 

/BARS/ 

Used  to  compute  0 for 

plate 

case , 

G 

BARG1 

/CONSTS/ 

Constant  portion  of 

G for 

given 

plate 

case 

BARH 

/BARS/ 

Used  to  compute  0 for 

plate 

case , 

H 

BARH1 

/CONSTS/ 

Constant  portion  of 

H for 

given 

plate 

case 

BARJ 

/BARS/ 

Used  to  compute  *<S  for 

plate 

case , 

blast 

load , J 

BARJ1 

/CONSTS/ 

Constant  portion  of 

J for 

given 

plate 

case 

BARK 

/BARS/ 

Used  to  compute  <5  for 

plate 

case , 

blast 

load , K 

BARR 

/BARS/ 

Used  to  compute  0 for 

plate 

case , 

linear 

load , R 

TABLE  A. 3. 1.1.  LIST  OF  PROGRAM  VARIABLES,  Continued 


’ 
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1 
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BARS 


BART 


BARTHE 


BARU 

BEAM 


*BETA 

BETACB 

BETAFR 

BETASQ 

BETAX 

BFTNX 


BLANK 


/BARS/  Used  to  compute  0 for  plate  case,  linear 
load , S 

/BARS/  Used  to  compute  6 for  plate  case,  linear 
load,  T 

/BARS/  Used  to  compute  0 for  plate  case: 

A*B,  if  blast  load 
R*S,  if  linear  load 

/BARS/  Used  to  compute  6 for  plate  case,  linear 
load,  U 

(GETITL)  Five  constant  characters,  "BEAM  ", 
used  to  match  first  five  columns  of  first  input 
cards  as  a Beam  case 

/INPUTS/  Spatial  Pressure  Decay  Constant,  B, 

dimensionless 

/CONSTS/  B3 

/CONSTS/  B4 

/CONSTS/  B2 

(COMP)  Used  to  consolidate  calculations,  B*CURRX 
(BFTNX ,CALXH0 ) Subfunction  and  value  of  function 
defining  Xq,  used  by  BISECT  in  determining  root 
for  Beam  case 

(TCNTRL)  Alphabetic  constant  containing  one 
blank  character,  used  to  initialize  the  Failure 


Flag  for  printout 
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TABLE  A. 3. 1.1.  LIST  OF  PROGRAM  VARIABLES,  Continued 


i 


BPFLAG 

BTHED2 

BTHED3 

BVDOT 

BWFDOT 

BWPDOT 

CURRX 

*D 

DELTA 

DELTAK 

DONEFG 


/FLAGS/  Beam-Plate-Flag  to  denote  type  of  case: 

1,  Beam  Case 

2,  Plate  Case 

/COMPS/  Constant  portion  of  second  part  ol  0 
computation  for  Beam  Case 

/COMPS/  Constant  portion  of  third  part  of  0 
computation  for  Beam  Case 

/COMPS/  Constant  portion  of  V(PE)  computation 
for  Beam  Case 

/COMPS/  Constant  portion  of  WF  computation  for 
Beam  Case 

/COMPS/  Constant  portion  of  WP  computation  for 
Beam  Case 

(COMP)  Value  of  X as  determined  by  current  Runge- 
Kutta  arguments  for  6 and  0. 

/INPUTS/  Distance  between  reinforcing  rod  and., 
edge  of  opposite  side,  d,  in. 

/RESULT/  Midpoint  deflection,  6,  in. 

/CONSTS/  Constant  used  by  several  equations, 
nw/m 

/FLAGS/  Case-is-Done-Flag  (integer):  Set  to 
one  when  a condition  is  found  which  indicates 
the  case  should  be  terminated  (including  a 
normal  end);  set  to  zero  otherwise 
94. 
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TABLE  A. 3. 1.1.  LIST  OF  PROGRAM  VARIABLES,  Continued 


DX  ( 5 ) 


DXDX ( 2 ) 


EOFLAG 


EPSILN 


EPSLNU 


ERRFLG 


EXPBCX 


EXPBET 


(RUNGEK)  Runge-Kutta  computed  first  differen- 
tial changes  for  given  time  step  for  the  five 
differential  equations 

(RUNGEK)  Runge-Kutta  computed  second  differen- 
tial changes  for  given  time  step  for  the  two 
second  order  differential  equations 
/FLAGS/  End-of-File-Flag  (integer):  Set  to 
one  when  input  is  expended;  zero  otherwise 
(BISECT)  Acceptable  accuracy  error  in  converg- 
ence used  by  root-finding  routine 
/CONSTS/  Ultimate  strain  of  reinforcing  element, 
assumed  to  be  a constant,  0.2,  eu,  dimensionless 
(BISECT)  Error-Flag  returned  to  calling  pro- 
cedure (integer): 

0 = no  error 

+1  = error:  even  number  of  roots  in  given 

interval,  could  not  use  method  to  find  root 

-1  = error:  method  would  not  converge  after 
100  tries,  root  may  be  acceptable 
(COMP)  expB(CURRX-l) , where  CURRX  is  the  current 
value  of  X being  used  in  the  Runge-Kutta  computation 
/CONSTS/  exp(-B) 
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TABLE  A. 3. 1.1.  LIST  OF  PROGRAM  VARIABLES,  Continued 

EXPBX1  /BARS/  expB(XX-l),  where  XX  is  current  value 

being  considered  for  X 

EXPBZX  /BARS/  expB(Z*XX-l)  , where  XX  is  current  value 

being  considered  for  X 

EXPBZ1  (FTNZ)  expS(TRIALZ-l) , where  TRIALZ  is  current 

guess  at  value  for  z 

*F  /INPUTS/  Support  factor: 

1 = simply  supported  edges 

2 = clamped  edges 

FLAG  /PRINTS/  Failure  Flag,  used  to  indicate  whether 

or  not  reinforcing  element  has  fractured  (alpha- 
betic): blank,  no  fracture;  *,  fracture 

FOURTH  /CONSTS/  Constant,  1/4 

FTN  (BISECT)  Dummy  parameter  naming  function  sub- 

procedure for  which  root  is  being  found 


r 


FTNT 

FTNZ 


FI 


(COMP)  f(t) 

(PTZERO,  FTNZ)  Subfunction  and  value  of  func- 
tion defining  z,  used  by  BISECT  in  determining 

root,  plate  case 

(BISECT)  Value  of  given  function  at  lower  end 
of  current  interval  of  guesses,  f(Xl) 

(BISECT)  Value  of  given  function  at  upper  end 
of  current  interval  of  guesses,  f(X2) 


y 6 . 


F2 


TABLE  A. 3. 1.1.  LIST  OF  PROGRAM  VARIABLES,  Continued 


F3 


GUESS1 


GUESS2 


*H 

HALF 

HALFTI 

I 


IERRFG 


ITER 

J 


KOUNT 


(BISECT)  Value  of  given  function  at  current  guess 
being  tested,  f(X3) 

(BISECT)  Dummy  variable  for  lower  end  of  inter- 
val containing  root  used  as  first  seed  for  Bi- 
section method  in  BISECT 

(BISECT)  Dummy  variable  for  upper  end  of  inter- 
val containing  root  used  as  second  seed  for  Bi- 
section method  in  BISECT 

/INPUTS/  Beam  or  Plate  thickness,  h,  in. 

/CONSTS/  Constant,  1/2 

(RUNGEK)  One-half  given  time  step,  dt/2 
(COMP,  PAGE,  RUNGEK)  Array  subscript  varying 
over  a loop 

/FLAGS/  Input-Error-Flag:  Set  to  one  when  bad 
data  is  encountered;  zero  otherwise 
(BISECT)  Current  interation 

(CONCRE)  Current  Computer  System  print  density 
setting  as  returned  from  System  Library  routine, 
PDEN,  lines/in. 

/FLAGS/  Current  total  for  run  of  accumulated 
errors  caused  by  improper  input  data,  used  by 
CDC  FORTRAN  Extended  Version  4 library  routine 
ERRSET , to  inhibit  job  termination  for  bad  data; 
maximum  arbitrarily  set  to  100  in  CONCRE 


TABLE  A. 3. 1.1.  LIST  OF  PROGRAM  VARIABLES,  Continued 

LCOUNT  (RWDATA)  Value  of  KOUNT  before  reading  data 

for  new  case 

LMBPSM  (COMP)  Coded  sum  of  the  Load,  Mechanism,  and 

Beam-Plate-Flags  plus  one,  used  to  determine 
which  of  eight  different  sets  of  differential 
equations  will  be  computed: 

1 Beam  Case  in  Mechanism  1,  Linear  Load 

2 Beam  Case  in  Mechanism  2,  Linear  Load 

3 Plate  Case  in  Mechanism  1,  Linear  Load 

4 Plate  Case  in  Mechanism  2 , Linear  Load 

5 Beam  Case  in  Mechanism  1,  Blast  Load 

6 Beam  Case  in  Mechanism  2,  Blast  Load 

7 Plate  Case  in  Mechanism  1,  Blast  Load 

8 Plate  Case  in  Mechanism  2,  Blast  Load 

LOADFG  /FLAGS/  Load-Flag: 

1,  if  case  in  under  a Linear  Load 

2 , if  case  in  under  a Blast  Load 
LOOPFG  (TCNTRL , TSTEP)  Loop-Control-Flag:  Set  to 

one  if  loop  for  time  = t^ (TINCR)  ,t 2 > . • • , 
tmax  (TMAX)  should  be  exited  for  any  reason, 
including  normal  termination;  zero  otherwise 
MAXLIN  /PRINTS/  Maximum  number  of  lines  :o  be  printed 

per  output  page;  determined  from  print  density 


W 


TABLE  A. 3. 1.1.  LIST  OF  PROGRAM  VARIABLES , Continued 


MAXTRI  (TCNTRL)  Maximum  number  of  times  the  loop  for 

time  = t, ,t  , ...,t  can  be  tried  using  success- 
j-  z max 

ively  smaller  time  steps,  if  exited  abnormally 
MECHFG  /FLAGS/  Mechanism-Flag: 

1 Case  in  Mechanism  1,  x^  = a 

2 Case  in  Mechanism  2 , 0 £ x^  < a 

M1STFG  /FLAGS/  First-Time-Flag  used  to  provide  a 

l 

dummy  value  of  CURRX  to  prime  the  first  two 

[passes  through  COMP  during  the  first  time  step: 

0 Have  completed  at  least  two  passes 

Kl 

through  COMP 

— 1 Have  completed  first  pass  through  COMP  t 

2 Have  not  yet  called  COMP 

NAIDA  /PRINTS/  Integer  value  of  n plus  two;  subscript 

for  array  of  type  of  slab  titles,  TYPSLB 
1 Horizontal  Slab  with  Explosive  above 
2 Vertical  Wall 

3  Horizontal  Slab  with  Explosive  below 
NCONFG  /FLAGS/  Non-Convergent-Root-Flag  used  to  signify 

if  root-finding  technique  found  a non-convergent 
or  out-of-bounds  root  for  the  given  function 
determining  x^  or  z: 


TABLE  A. 3.1.1.  LIST  OF  PROGRAM  VARIABLES,  Continued 


0 Root  found  is  valid 

1 Root  did  not  converge,  was  out-of-range, 
or  could  not  be  found  with  this  method 
using  seeds  (guesses)  given 

-1  Root  did  not  converge  after  100  inter- 
actions, but  may  be  close  enough  to  use 

NERR  (RWDATA)  Count  of  errors  found  by  reading  bad 

data  for  the  given  case;  i.e.,  input  errors 
accumulated  for  given  case 

NF  /PRINTS/  Integer  value  of  F,  Support  Factor; 

subscript  for  array  of  support  titles,  TYPSUP: 

1 Simply  supported  edges 

2 Clamped  edges 

NTRIES  (DRIVER,  TCNTRL , TZERO)  Number  of  times  loop 

for  computing  and  printing  output  lines  for 
t « tp  t2,..*»tmax  has  been  tried  for  given 
case  (<MAXTRI) 

NUMLIN  /PRINTS/  Number  of  lines  printed  on  current 

output  page  (<MAXLIN) 

NUMPAG  /PRINTS/  Page  number  of  current  output  page  for 

given  case 


100  . 


TABLE  A. 3.1.1.  LIST  OF  PROGRAM  VARIABLES,  Continued 


NWAVEF 


ONEMCX 

ONEMXX 

ONEMZ 

ONEMZX 

ONEPZ 

*PC 

*PE 

PFTNX 

PLATE 


/PRINTS/  Integer  value  of  WAVEFN,  Pressure 
Wave  Function;  used  as  subscript  for  Pressure 
Wave  Function  Title  array,  TYPWAV : 

1 General  time  function 

2 Square  Wave  Function 

(COMP)  1-CURRX,  where  CURRX  is  the  current  value 
of  X being  used  by  the  Runge-Kutta  computation 
(CALBAR)  1-XX,  where  XX  is  the  current  value  of 
X being  used  by  CALBAR 
/CONSTS/  1-z 

(CALBAR)  l-z*XX,  where  XX  is  the  current  value 
of  X being  used  by  CALBAR 
/CONSTS/  1+z 

/INPUTS/  Pressure  at  plate  or  bean  center  for 
varying  portion  of  the  load,  Pc,  psi 
/INPUTS/  Pressure  at  plate  or  beam  edges  for 
constant  load,  P^,  psi 

(CALXHO,  PFTNX)  Subfunction  and  value  of  func- 
tion defining  used  by  BISECT  in  determining 

root,  plate  case 

(GETITL)  Five  constant  characters,  "PLATE",  used 
to  match  first  five  columns  of  first  input  card 
as  a Plate  case 


101. 
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\ 

l 
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HI 


4 

■ 


I 
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TABLE  A. 3.1.1.  LIST  OF  PROGRAM  VARIABLES,  Continued 
SIXTH  /CONSTS/  Constant,  1/6 

STAR  (TSTEP)  Alphabetic  constant  containing  one 

asterisk,  used  to  set  the  Failure  Flag  for  print- 
out 


STEPCT 

T 

*TAU 

THETA 

THETAD 

THETAU 

THETU1 

THIRD 

TIMNOW 


"TINCR 
*TITLE ( 15 ) 


/PRINTS/  Count  of  number  of  time  steps  taken 
since  last  printout  line 

/RESULT/  Current  time  for  time  step,  t,  sec. 
/INPUTS/  Pressure  duration,  t,  sec. 

/RESULT/  Hinge  rotation,  6,  radians 
/RESULT/  First  differential  of  0 with  respect 
tot,  0 

(TSTEP)  Rotation  for  ultimate  strain  of  reinforc- 
ing element,  0U,  dimensionless 
/CONSTS/  Constant  portion  of  0U,  4heu/d^ 

/CONSTS/  Constant,  1/3 

/PRINTS/  Current  time  at  start  of  run  as  re- 
turned from  System  Library  routine,  TIME,  used 
in  printing  page  headings  to  define  run  output 
/INPUTS/  Time  step  increment,  dt , sec. 

/PRINTS/  Title  of  given  case  provided  from 
case  header  card  and  used  in  page  headings; 

15  x 5 characters  in  length 


1C  3 . 


-•  • 
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TABLE  A. 3. 1.1.  LIST  OF  PROGRAM  VARIABLES,  Continued 


*TMAX 


TODAY 


*TPRINT 


T PAT  10 


TRIALX 


TRIALZ 


TRXSQ 
TRZSQ 
*TYPE 
TYPLOAC  2 ) 


/INPUTS/  Maximum  allowable  value  for  t for 
given  case,  tmax,  sec. 

/PRINTS/  Current  date  at  start  of  run  as  re- 
turned from  System  Library  routine,  DATE:  used 
in  printing  page  headings  to  define  run  output 
/PRINTS/  Interval  of  time  allowed  between 
printed  time  steps,  tprint,  sec.:  tprint  £ dt , 
every  time  step  will  be  printed;  tpr^nt  > tm<?x, 
no  time  steps  will  be  printed;  dt  < tpriP^  <_  tmax, 
only  those  times  steps  which  are  rounded  multiples 
of  tprint  will  printed 
(COMP)  t/r,  dimensionless 

(BFTNX,  PFTNX)  Dummy  parameter  for  current  guess 
at  X 

(FTNZ)  Dummy  parameter  for  current  guess  at  z 
(BFTNX)  TRIALX2 
(FTNZ)  TRIALZ2 

/PRINTS/  Type  of  Case,  "BEAM  " or  "PLATE" 

(PAGE)  Array  containing  title  for  Type  of  Load, 
used  for  printing  page  headings;  subscripted  by 
LOADFG  (Load-Flag);  each  title  contains  1X6 
characters 
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TABLE A. 3.1.1.  LIST  OF  PROGRAM  VARIABLES,  Continued 


WF 

/RESULT/ 

Work  done  by  Pressure  Load,  WF , in-lbs 

WK 

/RESULT/ 

Kinetic  energy,  WK(KE),  in-lbs. 

WP 

/RESULT/ 

Plastic  energy,  WP,  in-- lbs. 

X 

/RESULT/ 

Ratio  of  current  hinge  location. 

V 

o 

fd 

x: 

X 

X < 1),  X,  dimensionless 

XH 

/RESULT/ 

Current  hinge  location  (0  < xg  < a), 

xh , in . 

XHO 

(CALXHO ) 

Original  hinge  location,  x^q , in. 

*XM 

/INPUTS/ 

2 3 

Mass  per  unit  area,  m,  Ibf-sec  /in 

XMU 

/CONSTS/ 

Hinge  moment  per  unit  length,  plate 

case,  Mu, 

lb-in/ in . 

, XX 

(CALBAR) 

Dummy  parameter  for  value  of  X used 

by  subroutine 

XXCUBE 

(CALBAR) 

XX3 

XX  SQ 

(CALBAR) 

XX2 

XO 

(CALXHO) 

Ratio  of  original  hinge  location, 

xh0/a’  xo 

, dimensionless 

XI 

(BISECT) 

Value  of  lower  end  of  interval  contain 

ing  current  guess  at  root 

X2 

(BISECT) 

Value  of  upper  end  of  interval  contain 

ing  current  guess  at  root 

X3 

■ 

(BISECT) 

Value  cf  midpoint  of  interval;  i.e., 

current  g 

uess  at  loot 

I 

, 


TABLE  A. 3.1.1.  LIST  OF  PROGRAM  VARIABLES,  Concluded 


YC  5 ) 


ZCUBE 


ZFOUR 


(RUNGEK)  Value  of  the  first  differentials,  6 
and  6,  and  the  variables,  WF , WP,  and  V(PE), 
after  the  Runge-Kutta  method  solves  their  first 
order  differential  equations: 

Y(l)  = 6 
Y(  2 ) = 0 
Y ( 3 ) = WF 
Y(  4 ) = WP 
Y( 5 ) = V ( PE ) 

/CONSTS/  Ratio  of  final  hinge  location  to  b, 

plate  case,  z,  dimensionless 

/CONSTS/  z3 

/CONSTS/  z4 

/CONSTS/  z2 


TABLE 

A . 3 . 1 . 2 . 

DIMENSIONS  OF  CONCRE  VARIABLE 

VARIABLE 

DIMENSION  (Contents) 

AA 

5 

(6,0,WF,WP,V(PE)) 

AO 

5 

( 6 , 0 j WF , WP , V ( PE ) ) 

A1 

5 

(6 ,0 ,WF,WP,V(P*E)) 

A2 

5 

(6,0,  WF , WP , V ( P*E ) ) 

A3 

5 

(6,0,  WF  , WP  , V ( P*E ) ) 

ARG 

5 

(dt  ,6 ,6 ,9,0) 

BX 

5 

( d6 , d0 , dWF , dWP , dV ( dPE ) ) 

DXDX 

2 

(d6  ,d0  ) 

TITLE 

15 

(15A5,  Case  Title) 

TYPLOA 

2 

(2A6,  Descriptive  Title) 

TYPSLB 

4, 

3 C 3(?4A8  ) , Descriptive  Title) 

TYPSUP 

2 

(2A7,  Descriptive  Title) 

TYPWAV 

2, 

2 (2(2A6),  Descriptive  Title) 

U 

2 

(6,0) 

Y 

5 

(6,0  ,VJF  ,WP  ,V  (PE) ) 

TABLE  A. 3. 2.1.  CCNCRE  PROCEDURES  AND  COMMON  BLOCKS 


X X X X X 


INPUTS 

(19) 

xxxxxx  xxxxxxxx 

CONSTS 

(30) 

X XXX  XXXX  XX 

COMPS 

(16) 

XX  X 

CO 

P4  f' 

X XXX 

< 1 1 

ca 

AjS 

XXXXXX XX XX  XX  XXXXX  X 

Pm 

,c 

4-> 

bO 

o o <u  / 
X O J/w 


w o!  to  j < j o o o h x.  k a: 

X W < H H K O Dh  ft!  2 CJ  X X P-.  bl  < X 

t ) > o h < H rr;  w w x u i 2 s w ui  a.  m x 

X H H H ,'j  X M [<1  N J W 2 tn  H H X X -J  UJ 

"5  « X i4  X O M [-<  H < M H X hvno  < X 

u Q X w a;  h H m a.  o ea  [m  cq  a,  E-<  aJ  o u c_> 


PAGE 


COMPS,  contains  constants  and  variables  used  by  the 
subroutine  COMP  to  compute  the  equations  of  work 


i 


■ 


r 


and  motion; 

CONSTS,  contains  miscellaneous  constants  and  vari- 
ables used  throughout  the  program; 

INPUTS,  contains  the  input  data  for  the  case; 

PRINTS,  contains  counts  and  alphabetic  information 
used  for  printing  the  output  headings,  etc;  and 
RESULT,  contains  the  numerical  results  computed  at 
each  time  step  for  printout. 

A further  description  of  the  content  and  definition  of  each 
COMMON  block  can  be  found  in  Table  A. 3.2.2. 

A. 3. 3 Core  Requirements 

The  CONCRE  program  consists  of  21  user  supplied  procedures 
using  approximately  1800  cards  (or  lines).  It  also  uses 
several  FORTRAN  and  System  Library  routines.  The  program 
compiles  in  less  than  65K  words  of  memory  and  can  be  run 
either  interactively  or  on  batch. 

A breakdown  of  the  core  required  for  each  procedure  is 
given  in  Table  A. 3. 3.1.  This  shows  a total  of  2080  core 
positions.  However,  using  labeled  COMMON  blocks  causes  the 
CDC  6600  N0S/BE1  Operating  System  to  set  aside  a minimum  of 
52K  memory  positions  for  COMMON;  so,  the  actual  core  required 
to  run  the  program  is  54K . 


TABLE  A. 3. 2. 2.  COMMON  BLOCK  DESCRIPTIONS 


J 


t 


Cornmon 

Block 

Label  Variable 

Length 

Descriptive 

Title 

Defined 

In 

BARS  BARA 

(all  real) 

1 

A 

CALBAR 

BARB 

1 

IT 

CALBAR 

BARC 

1 

( J 

CALBAR 

BARD 

1 

D 

CALBAR 

BARDEL 

1 

JK  or  TU 

CALBAR 

BARDNM 

1 

Denominator  of 
A ,D , H , R 

CALBAR 

BARG 

1 

G 

CALBAR 

BARH 

1 

H 

CALBAR 

BARJ 

1 

J 

CALBAR 

BARK 

1 

K 

CALBAR 

BARR 

1 

E 

CALBAR 

BARS 

1 

S 

CALBAR 

BART 

1 

T 

CALBAR 

BARTHE 

1 

AB  or  RS 

CALBAR 

BARU 

1 

U 

CALBAR 

EXPBX1 

1 

expB(X-l) 

CALBAR 

EXPBZ1 

1 

expB ( zX- 1 ) 

CALBAR 

COMPS  ATHEDB 

(all  real) 

1 

Constant  part 
of  0 

BTZERC 

ATHEDL 

1 

Constant  part 
of  0 

BTZEF.O 

ATHED1 

1 

Constant  part 
of  0 

FT ZERO 

.. ... . : . 


• 
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TABLE  A. 3. 2.2.  COMMON  BLOCK  DESCRIPTIONS,  Continued 


ATHED2 

1 

Constant 
of  0 

part 

BTZERO 

PTZERO 

AVDOT 

1 

Constant 
of  V ( P E ) 

part 

BTZERO 

PTZERO 

AWFDB 

1 

Constant 
of  WF 

part 

BTZERO 

PTZERO 

AV/FDL 

1 

Cons.tant 
of  WF 

part 

BTZERO 

PTZERO 

BTHED2 

1 

Constant 
of  0 

part 

BTZERO 

BTHED3 

1 

Constant 
of  0 

part 

BTZERO 

BVDOT 

1 

Constant 
of  VCPE) 

part 

BTZERO 

' 

BWFDOT 

1 

Cons.tant 
of  WF 

part 

BTZERO 

BWPDOT 

1 

Cons.tant 
of  WP 

part 

BTZERO 

PWFDB 

1 

Cons.tant 
of  WF 

part 

PTZERO 

PWFDL 

1 

Cons.tant 
of  WF 

part 

PTZERO 

PWPDOT 

1 

Cons.tant 
of  WP 

part 

PTZERO 

PVDOT 

1 

Con§ta/it 
of  VCPE) 

part 

PTZERO 

CONSTS 
(all  real) 

ACUBE 

1 

a3 

TZERO 

ARSQ 

1 

AR? 

TZERO 

TABLE  A. 3. 2. 2.  COMMON  BLOCK  DESCRIPTIONS,  Continued 


ARSZP1 

1 

( AR2 ) ( z ) + 1 

PTZERO 

ARZSP1 

1 

(AR) ( z2 ) +1 

PTZERO 

ASQ 

1 

a2 

TZERO 

E 

1 

b 

BTZERO 

PTZERO 

BARA1 

1 

Constant  part 
of  A 

PTZERO 

BAR  Cl 

1 

Constant  part 
of  C 

PTZERO 

BARG1 

1 

Constant  part 
of  G 

PTZERO 

BARH1 

1 

Constant  part 
of  H 

PTZERO 

BARJ1 

1 

Constant  part 
of  J 

PTZERO 

BETACB 

1 

B3 

TZERO 

BETAFR 

1 

34 

TZERO 

BETASQ 

1 

e2 

TZERO 

DELTAK 

1 

nw/m 

TZERO 

EPSLNU 

1 

0.2 

CONCRE 

EXPBET 

1 

exp(-3) 

TZERO 

FOURTH 

1 

1/4 

CONCRE 

HALF 

1 

1/2 

CONCRE 

ONEMZ 

1 

1-z 

PTZERO 

OUEPZ 

1 

1 + z 

PTZERO 

r 


TABLE  A. 3. 2. 2.  COMMON  BLOCK  DESCRIPTIONS,  Continued 


i 


SIXTH 

1 

1/6 

CONCRE 

THETU1 

1 

4heu/d  2 

TZERO 

THIRD 

1 

1/3 

CONCRE 

W 

1 

w=386 . 4m 

TZERO 

XMU 

1 

Mu 

TZERO 

Z 

n 

± 

z 

PTZERO 

ZCUBE 

1 

Z3 

PTZERO 

ZFO'JR 

z4 

PTZERO 

ZSQ 

l 

z2 

PTZERO 

FLAGS 

(all 

integer) 

BADXFG 

7 

- 

Bad-Hinge-Flag 

TCNTRL 

CHEKXH 

BPFLAG 

1 

Beam-or- Plate- 
Flag 

NXTCAS 

GETITL 

DONEFG 

1 

Case-is-Done- 

Flag 

DRIVER 

TCNTRL 

TSTEP, 

CHEKXH 

EOFLAG 

1 

End-of-File- 

Flag 

CONCRE 

GETITL 

RWDATA 

IERRFG 

*i 

Input-Error- 

Flag 

DRIVER 

RWDATA 

KOUNT 

jL 

System  Input 
Error  Count 

CONCRE 

LOADFG 

1 

Type-of- Load- 
Flag 

TCNTRL 

MECHFG 

1 

Mechanism- Flag 

TCNTRL, 

CHEKXH 

i 
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TABLE  A. 3. 2. 2.  COMMON  BLOCK  DESCRIPTIONS,  Continued 


r 


M 


I 


I 


;■ 

ii 


M1STFG 

1 

First-Tiine- 

thru-COMP-Flag 

TCNTRL 

COMP 

NCONFG 

1 

Non-Convergent- 

Flag 

TCNTRL 

PTZERO 

CALXHO 

INPUTS 
(all  real) 

A 

1 

a 

RWDATA 

AIDA 

1 

n 

RWDATA 

ALPHA 

1 

a 

RWDATA 

AR 

1 

AR 

RWDATA 

BETA 

1 

6 

RWDATA 

D 

1 

d 

RWDATA 

F 

1 

F 

RWDATA 

H 

1 

h 

RWDATA 

PC 

1 

pc 

RWDATA 

PE 

1 

PE 

RWDATA 

Q 

1 

q 

RWDATA 

SIGMAC 

1 

ac 

RWDATA 

SIGMAR 

1 

ar 

RWDATA 

TAU 

1 

T 

RWDATA 

TINCR 

1 

dt 

RWDATA 

TMAX 

1 

^max 

RWDATA 

TPRINT 

1 

tprint 

RWDATA 

XM 

1 

m 

RWDATA 

WAVEFN 

1 

115. 

Pressure  Wave 
Function  Code 

RWDATA 

TABLE  A. 3. 2. 2.  COMMON  BLOCK  DESCRIPTIONS,  Continued 


! 


I 


F 


PRINTS  FLAG 

1 

Failure-Flag 
(alphabetic , 

Al) 

TCNTRL 

TSTEP 

MAXLIN 

1 

Maximum  lines/ 
page  (integer) 

CONCRE 

NAIDA 

1 

n+2  (integer) 

RWDATA 

NF 

1 

F (integer) 

RWDATA 

NUMLIN 

1 

Current  line 
number  on  list- 
ing (integer) 

PRINTR 

PAGE 

NUMPAG 

1 

Current  page 
number  on  list- 
ing (integer) 

RWDATA 

PAGE 

NWAVEF 

1 

Pressure  Wave 
Function 
( integer) 

RWDATA 

STEPCT 

1 

Count  of  inte- 
gration steps 
since  last 
printed  line 
(integer) 

TSTEP, 

PRINTR 

TIMNOW 

1 

Time  of  run 
(alphabetic , 

A9  ) 

CONCRE 

TITLE 

15 

Case  title 
(alphabetic , 
15A5 ) 

GETITL 

TODAY 

1 

Date  of  run 
(alphabetic , 
A10) 

CONCRE 

TYPE 

1 

Type  of  Case 
(alrhabetic , 

A5 ) 

GETITL 

«_ 


1 
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TABLE  A. 3. 2.2.  COMMON  BLOCK  DESCRIPTIONS,  Concluded 


RESULT 
(all  real) 

DELTA 

1 

6 

TCNTRL 

TSTEP, 

RUNGEK 

T 

1 

Current  time 
step,  t 

TCNTRL 

TSTEP 

THETA 

1 

0 

TCNTRL 

RUNGEK 

THETAD 

1 

6 

TCNTRL 

RUNGEK 

V 

1 

Potential 
Engergy,  V(PE) 

TCNTRL 

RUNGEK 

VEL 

1 

Velocity,  6 

TCNTRL 

TSTEP, 

RUNGEK 

WF 

1 

WF 

TCNTRL 

RUNGEK 

WK 

1 

Kinetic  Energy, 
WK(KE) 

TCNTRL 

TSTEP 

WP 

1 

WP 

TCNTRL 

RUNGEK 

X 

1 

X 

TZERO , 
CHEKXH 

XH 

1 

Current  hinge 
location, 

TZERO, 

CHEKXH 

i 


TABLE  A. 3. 3.1.  CONCRE  CORE  REQUIREMENTS  BY  PROCEDURE 


i 


i 

:1 


i 


If' 


Procedure 

Name 

Program 

Length 

Buffer 

Length 

COMMON 

Length 

CONCRE 

192 

94  8 

66 

DRIVER 

12 

10 

NXTCAS 

2 

36 

GETITL 

29 

36 

RWDATA 

241 

55 

TCNTRL 

8 1 

66 

TZERO 

112 

70 

BTZERO 

5 

75 

PTZERO 

23 

92 

CALXHO 

52 

29 

BISECT 

20 

FTNZ 

19 

59 

BFTNX 

13 

59 

PFTNX 

7 

66 

TSTEP 

52 

96 

RUNGEK 

47 

40 

COMP 

19 

103 

CALBAR 

16 

76 

CHEKXH 

43 

40 

PRINTR 

19 

37 

PAGE 

109 

36 

TOTALS 

1040 

948 

129 

TABLE  A. 

3.3.2.  CONCRE  CORE  REQUIREMENTS 

Types  of  Cases 

Procedures  Eliminated 

Core  Required 

Beam  only 

PTZERO,  FTNZ,  PFTNX, 
CALBAR 

2052 

Plate  only 

BTZERO,  BFTNX 

2099 

Beam  and  Plate 

None 

2117 
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This  program  is  set  up  to  run  both  plate  and  beam  cases. 


If  one  type  is  eliminated,  some  of  the  FORTRAN  procedures 
may  also  be  eliminated.  Table  A. 3.3. 2 shows  this  information 
along  with  the  savings  in  core  requirements. 

A.  3.4  Library  Routines 

As  mentioned  earlier,  the  program  uses  several  System 
and  FORTRAN  Library  Routines.  Tables  A. 3. 4.1  list  these 
routines  with  a short  description.  It  includes  the  names, 
in  parenthesis,  of  the  variables  affected  by  these  routines, 
as  well  as  the  CONCRE  subprocedure ( s ) from  which  each  is 
called . 

A. 4 Input  and  Output 

Input  to  the  program  is  in  80-column  card  format, 
using  input  unit  number  5,  while  the  only  output  of  the 
program  is  a 133-character  listing  or  print  file,  using 
output  unit  number  6 . 

A. 4.1  Program  Input  Data 

The  six  ordered  input  card  formats  shown  in  Table 
A. 4. 1.1  describe  the  parameters  for  a single  case.  This 
table  includes  the  name,  in  parenthesis,  of  the  program 
variable  in  which  each  piece  of  input  data  is  stored. 
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TABLE  A. 3. 4.1.  DESCRIPTION  OF  SYSTEM  AND 


I 


J 

l 


FORTRAN  LIBRARY  ROUTINES  USED 

ABS(u)  provides  the  absolute  value  of  u;  a FORTRAN 

intrinsic  function  (EPSILN  in  BISECT) 

DATE(u)  returns  the  current  date  as  the  value  of  u 

in  alphabetic  form  (lOKbmm/dd/yyb ) ; a FORTRAN 
utility  subroutine  (TODAY  in  CONCRE) 

EOF(u)  tests  for  an  end-of-file  condition  on  file 

unit  u,  returning  the  real  value  0.0  if  not 
set  or  1.0  if  an  end-of-file  condition  was 
found;  a FORTRAN  utility  function  (file  5 in 
GETITL  and  RWDATA) 

ERRSET  (u-pi^)  sets  the  maximum  number  of  input  errors  to 
be  allowed  before  program  termination  to  U2 
and  keeps  the  count  of  such  errors  in  u^;  a 
FORTRAN  utility  subroutine  (KOUNT  in  CONCRE) 
EXP(u)  computes  eu ; a FORTRAN  external  function 

(EXPBET  in  TZERQ , EXPBZ1  in  FTNZ , BFTNZ  in 
BFTNZ , EXPBX1  and  EXPBZX  in  CALBAR,  and  FTNT 
and  EXPBCX  in  COMP) 

PDEN(i)  returns  the  current  print  density  setting 

(lines/inch)  for  the  printer  through  the  value 
of  the  argument  i;  a System;  Library  routine 
(J,  MAX LIN  in  CONCRE) 


TABLE  A. 3. 4.1. 


DESCRIPTION  OF  SYSTEM  AND 


FORTRAN  LIBRARY  ROUTINES  USED,  Concluded 


SIGN  (u1>U  ) 


SQRT  (u) 


TIME  (u) 


multiplies  the  sign  of  U2  by  the  value  of 
u-]_;  a FORTRAN  intrinsic  function  (SF1,  SF2, 

SF3  in  BISECT) 

takes  the  square  root  of  u;  a FORTRAN  external 
function  (THETAU  in  TSTEP) 

returns  the  current  clock  time  as  the  value 
of  the  argument  u in  alphabetic  form  ( 9Hbhh . mm. ss ) 
a FORTRAN  utility  subroutine  (TIMNOW  in  CONCRE) 
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TABLE  A. 4. 1.1.  DESCRIPTION  OF  INPUT  FOR  CONCRE 


Card  2: 


Card  3: 


Card  4: 


( 1 6 A 5 ) TYPE,  TITLE 

Type  of  case,  beginning  in  column  1, 
"BEAM  " or  "PLATE" 

Identifying  title  for  case,  free-fiel 

beginning  in  or  after  column  6 , up  to 

75  characters  allowed 

(4F12.1)  A,  H,  AR,  XM 

Flate  Half  Width  or  Beam  Half  Span, 

a , in. 

Thickness  of  Plate  or  Beam,  h,  in. 
Length  to  Width  Ratio,  AR,  dimension 

Mass  per  Unit  Area,  m,  lbs.  sec'V-in. 
(4F12.1)  PC,  PE,  ALPHA,  TAU 
Maximum  Distributed  Pressure  Load, 
Uniform  Pressure  Load,  P£ , psi 
Pressure  Decay  Constant,  a,  dimensions 
less 

Duration  of  Pressure,  t,  sec. 

(4 FI 2.1)  F,  WAVEFN , AIDA,  BETA 
Type  of  Support  for  Beam  or  Plate,  F, 
dimensionless : 

1,  Simply-Supported 

2,  Clamped-Supported 


(TYPE) 


(TITLE) 


(A) 

(H) 

ess  (A.R) 

(XM) 

psi  (PC) 
(PE) 

(ALPHA) 

(TAU) 

(F) 


TABLE  A. 4. 1.1. 


DESCRIPTION  Or  INPUT  TOR  CONCRE,  Concluded 


Card  5: 


Card  6: 


Type  of  Pressure  Wave  Function; 
dimensionless:  (WAVEFN) 

1,  General  Wave 

2,  Square  Wave  Function 

Weight  Vector,  n,  dimensionless:  (AIDA) 

0,  Vertical  Wall 

1,  Horizontal  Slab,  Explosive  Below 
-1,  Horizontal  Slab,  Explosive  Above 

Spatial  Pressure  Decay  Constant,  g, 
dimensionless  (BETA) 

(4F12.1)  SIGMAC,  SIGMAR,  Q,  D 
Concrete  Compressive  Strength,  oQ , psi  (SIGMAC) 
Reinforced  Steel  Yield  Stress,  or , psi  (SIGMAR) 
Reinforcement  Ratio  in  Tension,  q. 


dimensionless  (Q) 

Distance  Between  Tensile  Reinforcing  Rod 
and  opposite  edge,  d,  in.  (D) 

(3F12.1)  TINCR,  TMAX , TPRINT 

Integration  Time  Increment,  dt , sec.  (TINCR) 

Integration  Stop  Time,  tmax,  sec.  (TMAX) 

Integration  Time  Interval  between 
printed  lines,  tpran^. , sec.  (TPRINT) 
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Any  number  of  cases  can  be  run  one  after  the  other  by 
having  the  sets  of  six  cards  for  each  case  immediately 

follow  one  another.  Beam  and  Plate  cases  can  be  placed  in 
the  same  input  deck.  Errors  found  in  one  case  should  not 
affect  the  running  of  any  cases  which  follow  it.  (Exception 
if  a case  has  less  than  six  input  cards,  the  program  will 
attempt  to  read  some  of  the  cards  in  the  case  which  follows, 
and  will  skip  over  that  case. 

A. 4. 2 Program  Output 

Output  from  the  program  CONCRE  is  directed  to  the 
printer.  It  includes  listing  all  input  data  as  read,  the 
results  desired,  and  several  informative/error  messages. 

The  output  of  a sample  case  is  shown  in  Section  3.3. 

Table  A. 4. 2.1  shows  the  information  included  in  a 
normal  output  listing,  with  the  corresponding  program 
variable  names  given  in  parenthesis.  The  information  from 
the  first  input  card  (Type  and  Title)  is  included  in  all 
the  page  headings  for  the  case.  The  numerical  input  data 
is  printed  out  on  the  first  output  page,  as  well  as  some 
constants  computed  for  the  case.  The  second  and  succeed- 
ing pages  show  the  computed  results,  at  the  time  steps 
specified . 

If  errors  are  found  in  reading  or'  computing  initial 
case  variables,  only  the  first  page  is  printed  with  the 


TABLE  A. 4. 2.1.  CONCRE  OUTPUT 


» 


j 


Page  Headings  (printed  from  PAGE) 

Type  of  Case,  "BEAM  " or  "PLATE" 
Current  Date,  DD/MM/YY 
Time  at  start  of  run,  HH.MM.SS 
Case  Title 
Page  Number 

Type  of  Support  for  Slab 
Type  of  Pressure  Wave 

Position  of  Slab  relative  to  Explosive 
Type  of  Load 

First  Page  (Input  data  printed  from  RWDATA) 
a,  in. 
h,  in. 

Length  to  Width  Ratio,  dimensionless 
m,  lbf .-sec.  /in. 

Pc,  psi 
PE,  psi 

a,  dimensionless 
x , sec 

Support  Factor,  F 
Type  of  Pressure  Wave 
Weight  Vector,  n 


(TYPE) 
(TODAY) 
(TIMNOW) 
(TITLE) 
(NUMPAG ) 
(NF) 
(NWAVEF) 
(NAIDA) 
( LOADFG ) 

(A) 

(H) 

(AR) 

(XM) 

(PC) 

(PE) 

(ALPHA) 

(TAU) 

(F) 

( WAVEFN ) 
(AIDA) 


* 


TABLE  A. 4. 2.1.  CONCRE  OUTPUT,  Concluded 


S t 


B,  dimensionless  (BETA) 

crc , psi  (SIGMAC) 

ar,  psi  (SIGMAR) 

q,  dimensionless  (Q) 

d,  in.  (D) 

dt,  sec.  (TINCR) 

tmax  > sec>  (TMAX ) 

"tpj^Ln't5^ec.  ( T PRIN  T ) 

Second  and  Following  Pages  (Results  printed  from  PRINTR) 
Current  time  step,  t,  sec.  (T) 

Rotation,  6,  rad.  (THETA) 

Failure  Flag  (">•")  indicating  fracture  of 
reinforcing  element  (0>0U)  (FLAG) 

Midpoint  Velocity,  <5,  in. /sec.  (VEL) 

Midpoint  Deflection,  <5,  in.  (DELTA) 

Work  done  by  Pressure  Load,  WF , in. -lbs.  (WF) 

Plastic  Energy,  WP,  in. -lbs.  (WP) 

Kinetic  Energy,  WK,  in. -lbs.  (WK) 

Current  Hinge  Location,  x-.  , in.  (XH) 


12  f. 


appropriate  error  message 


many  error  messages  that  they  run  into  the  second  page 


but  no  headings  will  be  printed.)  Other  problems  can  be 


found  once  computations  are  underway;  these  will  be  printed 


terminated.  A description  of  all  the  error  message 


in  Table  A. 4. 2. 2.  This  include 


description  of  the  cause  of  the  message  with  possible  steps 


for  the  user  to  take,  and  the  name  of  the  FORTRAN 


cedure  in  parenthesis  from  which  the  message  was  printed 


A. 5 Numerical  Techniques 


ection  describe 


the  numerical  techniques  use 


by  the  program,  CONCRE.  These  are  the  bisection  (or  binary 


search)  method  used  by  the  subroutine  BISECT  to  find  the 


roots  of  the  equations  defining  Xq  and  z and  the  fourth- 
order  Runge-Kutta  technique  used  to  solve  the  five  simul 


taneous  differential  equations  of  work  and  motion  in  the 


RUNGE-KUTTA  function 


A. 5.1  The  Bisection  Method 


The  bisection  (or  binary  search)  method  uses  the  fact 
that  a function  which  crosses  the  x-axis  (or  goes  to  zero) 


ithin  a given  interval  must  have  different  signs  at  the  end 


points  of  that  interval.  Because  of  this  assumption,  thi 


TABLE  A.4.2.2.  CONCRE  ERROR  MESSAGES 

1.  INPUT  ERRORS  FOUND- -CASE  TERMINATED. 

Invalid  or  non-numeric  data  has  been  encountered  while 
reading  numeric  data  for  this  case.  This  message  will 
follow  other  messages  which  should  better  specify  the 
error(s).  The  program  will  proceed  to  the  next  case. 
(DRIVER) 

2.  BAD  TYPE  £ TITLE  CARD  FOUND 

xxxxxxxxxxxxxxxxx  . . . xxxxxx 

CARDS  MAY  BE  OUT  OF  ORDER 

The  program  was  trying  to  read  the  Type  £ Title  (or 
first)  card  of  a new  case,  but  found  a car'd  which  did 

not  contain  "BEAM  " or  "PLATE"  in  the  first  five  columns. 
The  card  found  is  printed  in  the  message  and  is  consider- 
ed erroneous.  The  program  will  try  to  read  the  next  card 
as  a Type  £ Title  card.  (GETITL) 

3.  END  OF  FILE  FOUND  DURING  INPUT  OF  DATA 

CASE  IS  TERMINATED  AS  NOT  ALL  DATA  IS  PRESENT 
An  end-of-file  condition  was  encountered  while  the  pro- 
gram was  trying  to  read  the  num.eric  data  (cards  2 through 
G)  for  the  giver  case.  The  case  is  terminated  and  the 
program  will  • rint  Message  5 and  begin  normal  termination. 
( RWDATA ) 


TABLE  A. 4. 2. 2.  CONCRE  ERROR  MESSAGES,  Continued 


4.  DATA  ERROR ( S ) WERE  FOUND  IN  INPUT 

CASE  IS  TERMINATED  AS  DATA  IS  IN  INCORRECT  FORM 
Non-numeric  data  has  been  encountered  while  reading  in- 
put cards  2 through  6 for  the  given  case.  This  message 
will  follow  System  Error  messages  which  should  show  the 
erroneous  card(s).  The  program  will  print  Message  5 and 
Message  1 and  proceed  to  the  next  case.  (RWDATA) 

5.  SOME  OF  THE  VALUES  LISTED  ABOVE  MAY  ACTUALLY  BE  FROM 

PREVIOUS  CASE 

The  input  data  for  the  case  has  been  printed  on  this 
same  page.  However,  if  invalid  data  was  read,  the  System 
did  not  allow  the  new  invalid  value  for  the  given  var- 
iable to  replace  the  value  used  by  the  last  case.  This 
message  follows  Message  3 or  Message  4.  (RWDATA) 

6.  VALUES  OF  AIDA,  F,  OR  WAVEFN  ARE  OUT  OF  RANGE 
CASE  IS  TERMINATED. 

The  correct  values  are  1,  0,  or  -1  for  AIDA  (the  Weight 
Vector),  1 or  2 for  F (the  Support  Factor),  and  1 or  2 
for  WAVEFN  (the  Wave  Function).  One  or  more  of  these  is 
not  correct.  The  program  will  print  Message  1,  and  begin 
on  a new  case.  (RWDATA) 

7.  SOME  COMPUTED  CONSTANTS  FOR  THIS  CASE  ARE  NEGATIVE 
SOME  INPUT  VALUES  MUST  BE  IN  ERROR 
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TABLE  A. 4. 2. 2.  CONCRE  ERROR  MESSAGES,  Continued 


CASE  IS  TERMINATED 

The  values  of  B (Beam  Width  or  Plate  Half  Length,  b), 

XMU  (Hinge  Moment,  Mu),  and  W (Weight  per  unit  Area,  w) 
are  computed  from  input  data  values  and  should  not  be 
negative.  Recheck  input  values.  The  program  will  pro- 
ceed to  the  next  case.  (TZERO) 

8.  THE  ROOT  OF  THE  EQUATION  FOR  Z WAS  NOT  BETWEEN  ZERO  AND 

ONE 

The  root-finding  procedure,  BISECT,  failed  to  find  a root 
for  z (the  Ratio  of  the  Final  Hinge  Location  to  the  Plate 
Half  Length)  in  the  interval  (0,  1).  Some  of  the  input 
data  must  be  incorrect.  The  program  proceeds  to  the  next 
case.  ( PTZERO ) 

9.  NO  VALUE  OF  THE  ORIGINAL  HINGE  LOCATION,  XHO , WAS  FOUND 

IN  THE  INTERVAL  (0,  A). 

IT  IS  ASSUMED  TO  BE  THE  VALUE  OF  A. 

The  root-finding  procedure,  BISECT,  was  unable  to  find  a 
value  of  x^q  in  the  interval  (0,  A).  Our  experience  shows 
it  is  acceptable  to  assume  the  value  of  x ^ as  a.  (CALXHO) 

10.  THE  BISECTION  METHOD  USED  TO  FIND  THE  ORIGINAL  HINGE  LO- 
CATION, XHO,  DID  NOT  CONVERGE  AFTER  100  ITERATIONS. 

THE  RESULT  OF  THE  LAST  ITERATION  WILL  BE  ASSUMED  CORRECT: 
XHO  = XXX. 
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TABLE  A. 4. 2. 2.  CONCRE  ERROR  MESSAGES,  Continued 


■ 


J 


i 


The  root-finding  procedure,  BISECT,  failed  to  converge 
within  the  given  e (-0.00000001)  after  100  iterations. 
While  this  could  be  caused  by  a discontinuity  of  the 
given  function,  no  known  discontinuities  exist  in  the 
given  interval  for  either  the  Beam  or  Plate  functions. 

The  program  assumes  convergence  is  close,  sets  Xq  to  the 
average  of  the  last  two  values  used  by  the  Bisection 
method,  and  continues.  (CALXH0) 

11.  TIME  EXCEEDED 

The  maximum  integration  time  (TMAX)  has  been  reached 
without  finding  the  Maximum  Deflection  for  the  case. 

The  program  goes  to  the  next  case.  Increase  the  inte- 
gration stop  time  (TMAX)  and  rerun  case  to  obtain  fur- 
ther results.  (TSTEP) 

12.  INSUFFICIENT  PRESSURE  TO  GIVE  A RESPONSE 

The  Maximum  Deflection  was  reached  before  the  first 
integration  time  step  was  completed.  The  program  goes 
to  the  next  case.  If  input  data  looks  correct,  decrease 
value  of  the  integration  time  step  (TINCR)  and  rerun  the 
case.  (TSTEP) 

13.  HINGE  LOCATION  HAS  OVERSHOT  FINAL  HINGE  LOC 

FINAL  HINGE  LOC  = XXX  HINGE  IS  AT  XXX 


TAELE  A.  ‘1.2. 2. 


COHCRE  ERROR  MESSAGES,  Concluded 


i 
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TIME  INCREMENT  HAS  BEEN  HALVED--CASE  WILL  BE  RERUN 
Between  the  last  two  integration  time  steps,  the  hinge 
location,  x^,  moved  from  a Mechanism  2 position  past  the 
final  hinge  location  (more  than  2%  pact).  This  infers 
that  the  computations  are  too  inexact.  Hence,  the  program 
automatically  divides  the  integration  time  step  (TINCR)  by 
two  and  reruns  the  case,  starting  again  at  t=0.  This  will 
only  be  done  a maximum  of  once  (MAXTRI  = 2);  then  the  case 
will  be  terminated.  (CHEKXH) 

14.  HINGE  LOCATION  IS  NEGATIVE- -CASE  IS  TERMINATED 
CHECK  INPUT  VALUES 

The  results  of  the  last  integration  step  have  produced 
a negative  hinge  location  (x^  = 6/0).  This  makes  no  sense; 
data  values  must  be  wrong.  (CHEKXH) 

15.  AN  ASTERISK  INDICATES'  THAT  A REINFORCING  ELEMENT  HAS 

FRACTURED 

This  is  a standard  footnote  printed  on  the  second  and 
succeeding  pages  of  the  output  for  every  case.  Should 
there  be  an  asterisk  following  the  value  printed  for  0, 
then  0 is  greater  than  0U  which  infers  the  failure  of  a 
reinforcing  rod.  (TSTEP,  PAGE) 


I 
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method  will  not  work  ior  finding  roots  of  functions  which 
only  become  tangent  to  the  x-axis.  Nor  will  it  work  for  a 
function  with  an  even  number  of  roots  within  the  given  inter- 
val, since  the  signs  of  the  function  at  the  ends  of  the  inter- 
val would  be  the  same. 

Basically,  the  method  works  in  the  following  manner. 
First,  it  checks  the  signs  of  the  function  at  the  end  points. 
If  they  are  the  same,  the  routine  terminates  in  error.  If 
they  are  not  the  same,  the  midpoint  of  the  interval  is  chosen 
and  the  sign  and  value  of  the  function  at  that  point  are  de- 
termined. If  the  value  is  zero,  the  root  has  been  found.  If 
not,  the  interval  is  halved  using  the  midpoint  to  replace 
the  previous  end  point  whose  function  had  the  same  sign. 

The  midpoint  of  the  new  interval  is  chosen  and  the  method 
is  repeated  until  the  function  or  the  interval  goes  to  zero. 

In  other  words,  the  interval  containing  the  root  is  made 
smaller  and  smaller  until  the  root  is  the  midpoint  of  that 
interval . 

The  three  functions  which  require  roots  found  by  this 
method  within  the  program  are  the  beam  equation  for  Xq  in 
the  interval  (0,  1)  (the  function  BFTNX),  the  plate  equation 
for  Xq  in  the  interval  (0,1)  (the  function  PFTNX),  and  the 
equation  for  z used  by  plate  cases  in  the  interval  [0,  1] 

(the  function  FTNZ).  Tests  run  during  the  formation  of  the 


program  showed  that  each  of  these  functions  had  one  real 
root,  at  most,  in  the  intervals  given. 

Many  other  root-finding  techniques  could  have  been 
used  instead  but  the  bisection  method  was  chosen  for  its 
simplicity  and  rapid  convergence. 


A. 5. 2 The  Fourth-Order  Runge-Kutta  Technique 

The  fourth-order  Runge-Kutta  technique  is  used  to 
solve  simultaneously  at  each  time  step  the  two  second- 
order  differential  equations  of  motion  (0  and  6)  with  the 

three  first  order  differential  equations  describing  the  work 

• • • * 

functions  (WF,  WP,  and  V(PL')).  These  five  equations  are 
functions  of  time,  t,  the  hinge  location,  x^ , and  the  veloc- 
ity of  rotation,  0 as  follows: 

(fp(t,x^),  Mechanism  2 
C , Mechanism  1 
0 = f2(t,xh) 

WF  = f 3 ( t , xh , 0 ) 

WP  = f 4 ( 0 ) 

V ( PE ) = f5(xh,6) 

Now  the  hinge  location,  x^ , is  defined  in  terms  of  0 and  6 
as : 

xh  = ^ / 9 

So,  our  equations  become,  in  Mechanism  2 (the  more  complex 


form)  : 


6 = f -j  (t  ,6 ,6) 

0 = f2(t,fi,0) 

WF  = f3(t,6,e,0) 

WP  = f4(0) 

v ( pe ) = f5(fi,e,e) 

For  the  sa> e of  uniformity,  the  subprocedures  RUNGEK  and 
COMP  expect  each  of  these  five  functions  to  have  five  argu- 
ments; i.e.,  f p (t ,6 , 6 , 9 , 9 ) , even  though  6 is  never  used. 
When  the  Runge-Kutta  computations  for  the  time  step  have 
been  completed,  the  new  value  of  is  computed  from  the 
just-obtained  values  of  0 and  6 (in  CHEKXH) . 

Hence,  given  the  initial  conditions  and  letting  h be 
the  time  increment,  dt  (or  TIMCR)  , we  can  determine,  induc- 
tively, using  the  fourth-order  Runge-Kutta  technique  for 
simultaneous  second-order  equations, 

6t  + ] = <5-^ +h6t-Kh/6  ) (aQ+a2  + a2  ) 

6^.  + 2 = <5 ^ + ( 1/ 6 ) ( ag  + 2aq+ 2a2+S3  ) 

et+1  = 0t+h0t+(h/6) (b0+b!+b2) 

®t+l  = $t+ ^ 1/6 ) (bQ+2bi+2b2+b3 ) 

where 

a0  = hf 2 (t, 6, 6, 0,0) 

al  = hf  ^ (t+hh , 6+\h5 , 6+%3q  , G+^hO  , G+^bg ) 

a2  = hf]_  (t+lih,  6+1-ih6+’<haQ  , 6+^a-,  , O+^hG+^hbQ  , 0+^b^  ) 

a3  = hf 2 (t+h , 6+h6+hha2 , 6+a2 , G+hG+^hb^ , 0+b^ ) 
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bg  = hf ^ ( t , 6 , 6 , G , 0 ) 

= hf  ^ ( t+'ih  , 6 + ’;h6  , G+’ih G , G+]'ib^  ) 

bj  = hf  2 ( t + 3sh , G+JghS+^hag  , 6+ha^  , G+^hG  + ^hbg  , 0 + 3sb^ ) 

• • • • 
b3  = hf  2 (t+h  , 6+h6+%ha^  , <5  + ^2  5 G+b0+%hb-^ , 0 + b2  ) 

The  first-order  differential  equations  are  computed  at  the 

same  time ; e . g. , 

WFt+i  = WFt+(l/6)(cQ+2c1+2c2+c3) 

where 

Cq  = hf 3 ( t , 6 , 6 , 0 , 0 ) 

• • • • 

Cj  = hf  3 ( t+^h , 6+JghS  , (S+^ag  , 8+^h0 , 0+^bg  ) 

C2  = hf  3 ( t+^h  , 6 + '5h6+^hag  , 6+^a^  , 0+%h0+%hbg  , G+^b^  ) 

C3  = hf  3 ( t+h  , 6+h6+^haq  , 6 +a2 , 0 + hO+^hb-^  , 0 + b2  ) 
and  where  ag , aq,  a2,  bg,  bp,  and  b2  are  defined  as  above. 

The  equations  for  WP  and  V(PE)  are  similar. 

From  the  above,  it  is  clear  that  the  arguments  for  the 
five  different  equations  are  the  same  at  each  step,  i=0,l,2,3. 
These  are  computed  in  the  subprocedure  RUNGEK  and  passed  in 
the  array  ARG  to  the  subprocedure  COMP  which  actually  uses 
these  arguments  to  compute  ap , bp,  etc  for  the  given  step,  i. 
These  results  are  stored  in  the  arrays  AO,  Al,  A2,  and  A3 
through  the  dummy  parameter  array  AA. 

The  initial  conditions  for  each  differential  equation 
is  determined  earlier  in  the  program;  they  are  all  initialized 
to  zero  in  the  procedure  TCNTRL.  Since  the  hinge  location, 
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is  the  result  of  a division  of  0 which  would  be  zero  at  time 
zero,  we  alter  the  arguments  in  the  array  ARG  given  above 


Runge-Kutta  time  step;  fn(  t ,xhg  , <S , 1 . 0 , 0 ) . This  can  be  done 
as  the  two  arguments  affected,  <5  and  8,  are  used  solely  to 


compute  x^.  The  arguments  for  6 for  the  third  and  final 
parts  of  the  first  time  step  use  the  results  fnom  the  first 
two  parts  of  the  first  time  step  (stored  in  AO  and  Al)  and 
will  be  non-zero,  so  that  x^  can  be  determined  properly. 

An  alternate  approach  would  have  been  to  define  each 
equation  in  three  unknowns:  fn(t,x^>9).  In  this  notation 
the  equation  for  <$-£+1  would  contain  ag  through  ag  defined 
as  follows: 

ag  = hf^(  t ,Xj1  = 6/0 , 0 ) 

• • • 

a-^  = hf  j ( t+%h , ( 6+5jh6  ) / ( 6+^h0  ) , 0+5gbg  ) 

• • • 

a 2 = hf  2 ( t+^h  , ( 6+5gh6+lchaQ  ) / ( 0 + 'ih0+khbg  ) , O+^b^  ) 

a3  = h-fp  j (t+h  , ( 6+h6+%ha^  ) / ( B+hQ+hhb^  ) , 0+b2  ) , 
for  the  same  initial  conditions.  Again,  care  would  have  to 
be  taken  with  x^ , during  the  first  two  parts  of  the  first 


time  step. 
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A. 6.  PROGRAM  LISTING 


PkOGRA  M CGNCRE ( I NF  UT , OUTPUT , f APE5=iNPUT , T AP  ES  = U U f PUT ) 

C 

C WRITTEN  BY  C.  C.  SCHAU3LE  FOR  C.  A.  ROSS 

C UNIVERSITY  OF  FLORIDA  GFAOUATE  ENGINEERING  CENTER 

C EGLIN  A 1 K FORCE  BASE,  FLORIDA 

C MARCH,  1979 

C 

C THIS  PROGRAM,  CONCRE,  COMPUTES  THE  DEFLECTION  OF  THE  CENTER 

C OF  A REINFORCED  CONCRETE  BEAM  OR  FlATE  CAUSED  BY  A LINEAR 

C OR  BLAST  LOAD. 

C 

C THIS  MAIN  PROCEDURE  CALLS  THE  SYSTEM  SUBPROCEDURES,  DATE, 

C TIME,  PDEN,  AND  ERRSET  TO  INITIALIZE  THE  CURRENT  DATE, 

C TIME  AND  PRINT  DENSITY  AND  TO  OVERRIDE  THE  AUTOMATIC 

C TERMINATION  CAUSED  BY  INVALID  INPUT  DATA.  IT  ALSO  CALLS 

C THE  SUBPROCEDURE,  DRIVER,  TO  HANDLE  THE  READING  AND 

C PROCESSING  OF  EACH  CASE. 

C 

C TAPE5  IS  THE  NORMAL  INPUT  FILE  IN  8G-COLJMN  CARD  FORMAT. 


C TAPE6  IS  THE  LISTING  OR  PRINT  FILE. 

C 
C 

INTEGER  BAOXFG, 

COMMON  /FLAGS  / 

1 

COMMON  /PRINTS/ 

l 1 

COMMON  /CONSTS/ 

1 

2 

3 

4 
C 

C INITIALIZATION 

P 

Vy 

C SET  DATE  AND  TIME 

C GET  MAXIMUM  NUMBER  OF  LINES  FER  FAGE 

C SET  CONTROL  TO  HANDLE  INPUT  ERRORS 

C CLEAR  END-OF-FILE  FLAG 

C 


3PFLAG,  DCNEFG , EOFLAG 

BAOXFG,  BFFLAG,  DONEFG,  EOFLAG,  IERRFG,  KOUNT, 
LOAOFG,  MECHFG,  M1STFG,  NCCNFG 

FlAG,  MAXLIN,  NAIUA,  NF,  NUMLIN,  NUMPAG , NWAVEF, 
STEPCT , TIMNOW,  T 1 T L E ( 1 5) , TODAY,  TYPE 
ACUdE,  ARSQ ♦ ARSZP1,  ARZSF1,  AS3,  3 , OARA1, 
BAKC1,  BAFG1 , BARH1,  BARJ1 , 3ETACB,  BETAFR, 

BET  AS  Q,  DELTAK,  EPSlNU,  EXPBET,  FOURTH,  HALF, 
GNEMZ , CNEPZ,  SIXTH,  THETU1,  THIRD,  W,  XMU,  Z, 

Z CUBE  , ZFCUR,  ZSQ 


OOOO  O OOO  OOOOO  OO  OOOO  OOOOOOOO  ooo 


c 

c 

c 

r* 

L» 

C 

C 

C 

C 

C 

C 


SUB)  UJT  I NE  DRIVER 

HIS  SU3R0U1INE  DRIVES  ON  E CASE  FROM  BEGINNING  (READING 
0-  DATA)  TO  END  (COMPLETION  OF  EXECUTION  OK  ERRONEOUS 
TERMINATION)  . II  C ALi-S  SUBROUTINES  NxTCAS,  10  HANDLE  THE 
READING  CF  THE  DATA  FOR  THE  NEXT  CASE,  AND  TCNTRL,  TO 
CONTROL  THE  CASE  THROUGH  THE  TIME  STEP  PROCESSING,  FROM 
INITIALIZING  FOR  TIME  ZERO  THROUGH  PROCESSING  EACH  TIME 
STEP  UNTIL  TIME-MAX  OR  UNTIL  NORMAL  COMPLETION. 

H IS  CALLED  3 V THE  MAIN  PROGRAM,  CONCRE. 

INTEGER  BAJXFG,  BPFLAG,  DCNEFG,  EuFLAG 

COMMON  /FLAGS  / 6A0XFG,  BFFlAG,  DGNEFG,  EOFLAG,  IERRFG,  KOUNT , 
1 LOADFG,  KECHFG,  MIST  FG,  NCCNrG 

clear  input -error-flag  and  read  next  case 

IEKkFG  = 0 
CAlL  NXTCAS 


CHECK  INPUT  OR  END-OF-FILE  ERRORS 
IF  END  OF  FIlE  FOUND 
THEN  TERMINATE  CASE 
IF  - (EOFLAG. NE. G)  GO  TO  2000 
ELSE 

if  input  error  found 

IF  ( IERP.FG.EQ.  0)  GO  TO  5G0 

THEN  PRINT  MESSAGE  AND  PROCEED  TO  TERMINATE  CASE 
WRITE  (6,9901) 

GO  TO  1999 

ELSE  CONTINUE 

CLEAR  CASE-IS-DONE-PLAG  AND  NUM3ER-Or-TKI ES-COUNT 

5 0 C CONEFG  = 0 
NTrIES  = 0 

DO  WHILE  CASE-IS -CONE-FLAG  IS  CLEAR 

10G0  IF  (.NOT.  (DONEFG.EQ.G)  ) GO  TO  1999 

CALL  TCNTRL  (NTRIES) 

GJ  TO  100G 

END  WHILE  (GONEFG) 

END  IF  ( i E k R f G ) 

140. 


- .... 


1999  GO  TO  2000 


E >1 0 IF  (EOFLAG) 

2000  RLTURM 

FORMAT  STATEMENT 

9901  FORMAT  <lHOt*INPUT  ERRORS  FOLNC--CASE  T ERMi N AT EJ* ) 


ooo  ooooo  o oooo  oon 


subroutine  nxtcas 


C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


THIS  SUBROUTINE  OBTAINS  THE  INPUT  FOR  THE  N EX  7 CASE  AMD 
SETS  UP  THE  OUTPUT  HEADINGS  AND  LISTS  THE  INPUT. 

THE  OEAM-PL ATE-FLAG  (6PFLAG)  IS  SET  TO  SHOW  WHETHER  THE  CASE 
IS  A BEAM  OR  A PLATE  CASE.  THE  ENO-GF-FILE  FLAG  (EOFLAG) 

IS  SET  IF  AN  END-OF-FILE  INDICATOR  IS  FOUND  WHEN  DATA  IS 
EXPECTED. 

THIS  ROUTINE  CALLS  INC  SUBROUTINES! 

GETITL!  TO  LOOK  FOR  THE  FIRST  CARD  OF  INPUT  FOR  THE  CASE, 
A TYPUTITLE  CARD  (ALPHABETIC  INPUT) 

RWOATA!  TO  READ  AM  PRINT  OUT  THE  NUMERIC  INPUT  VALUES, 

AS  WELL  AS  ANY  INPUT  ERRORS. 

THIS  SUBROUTINE  IS  CALLED  BY  THE  ROUTINE,  DRIVER. 


INlEoEk  BAQXFG » 
COMMON  /FLAGS  / 

1 

COMMON  /PRINTS/ 

1 


GFFLAG,  CONEFG,  EOFLAG 

BADXFG,  BPFLAG,  OCNEFG , EOFLAG,  IERP.FG,  KOUNT, 
LOADFG,  M E 0 r F G , MlSTFG , NCONFG 

FLAG,  MAXLIN,  NAIjA,  NF,  NUMLIN,  NUPPAG,  NWAVE-, 
STFPCT,  TIMNOW,  TITLE(ID),  TODAY,  TYPE 


CLEAR  THE  BEAM- OR -PLATE  FLAG 


BFFLA3  = 0 

do  while  the  end-of-file  and  beam-or-pl ate  flags  are  clear 

(I.E.,  UNTIL  THE  FIRST  CARD  CF  a new  CASE  IS  READ) 

IOC  IF  (.NOT.  (DFFLAG.EQ.O  .AND.  EOFLAG. EQ.O)  ) 

1 GO  TO  1GC0 


CAL.  GETITL 

GO  TO  100 

END  WHILE  ( BPF LAG  AND  EOFLAG) 

IF  THE  END-OF-FILE  HAS  NOT  YET  BEEN  REACHED 

1000  IF  (EOFLAG. NE. 0)  GO  TO  5000 

THEN  READ  AND  FkINT  DATA  FOR  THE  NEW  CASE 

CALu  RwDA  T A 
C 

C ELSE  CONTINUE 

C END  IF  (NOT  ELI) 

C 

5 0 0 i RETURN 
C 

END 


SUbkOJTINE  GElAlL 


I 

i 


j 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

o 

Hi  0 c 
c 
c 


c 

c 

2 0 0 C 

c 


THIS  SUBROUTINE  SEARCHES  THE  INPUT  FILE  FOR  THE  NEXT 
CARD  (OR  RFCOkO)  WHICH  CONTAINS  THE  CHARACTERS  "BEAM  ” 

OR  ••PLATE**  IN  THE  FIRST  FIVE  POSITIONS  Or  THE  CARD, 
INDICATING  A TYPE  & TITLE  INFUT  CARD,  THE  FIRST  INPUT 
CARO  FOR  A CASE.  WHEN  SUCH  A CARD  IS  FOUND,  THE  INPUT 
INFORMATION  IS  STORED  IN  THE  COMMON  BLOCK,  PRINTS,  TO  BE 
USED  LATER  IN  PRINTING  OUTPUT  HEADINGS  FOR  THE  CASE. 

ANY  CARDS  FOUND  WITHOUT  "BEAM  " OR  "PLATE"  IN  THE  FIRST 
FIVE  COLUMNS  ARE  CONSIDERED  FRkONFOUS  DR  OUT-OF-ORDER  AND 
ARE  PRINTED  OUT  AS  ERRORS. 

THE  END-OF-FILE  FLAG  IS  SET  IF  AN  ENQ-OF-FILE  INDICATOR 
IS  DETECTED.  THE  BE  AM -OR- PL  AT E - F L AG  IS  SET  ACCORDING  TO  THE 
TYPE  uF  TYFEiTITLE  CARD  FOUND. 

THIS  SUBROUTINE  IS  CALLED  BY  THE  PROCEDURE  NXTCAS.  IT 
USES  THE  SYSTEM  SUBPROCEOURE , EOF,  TO  CHECK  FOR  AN 
END-OF-FILE  CONDITION. 


INTEGER  BADXFG,  BPFLAG,  DON  EFG , EOFLAG 

COMMON  /FLAGS  / BADXFG,  BPFLAG,  DON  EF  G,  EOF. AG,  1ERRFG,  KOUNT, 

1 LOADFG,  MECHFG,  MlSTFG,  NCONFG 

COMMON  /PRINTS/  FLAG,  MAXLIN,  N A I U A , NF , NUMLIN,  NUMPAG,  NWAVEF, 
1 STEPCT , TIMNOW,  T1TLE(15),  TODAY,  TYPE 

DATA  BEAM  /5H8EAM  /, 

1 PLATE  /5HPLATE/ 


READ  NEXT  RECORO 

REAO  (5,9900)  TYPE,  TITLE 

IF  END-OF-FILE  HAS  BEEN  REACHED 

IF  (EOF (5) .NE. 1.0)  GO  TO  1000 
THEN  SET  END-OF-FILE  FLAG 
EOF. AG  = 1 
GO  TO  9000 

ELSE  CONTINUE  BY  CHECKING  TYPE 
CASE  ON  TYPE 

IF  THIS  IS  A BEAM  CASE 
IF  (TYPE.NE.bEAM)  GO  TO  2000 

THEN  SET  bEAM-PLATE  FLAG  TO  INDICATE  BEAM  (1) 
BPFLAG  = 1 
GO  TO  9000 

IF  THIS  IS  A PLATE  CASE 
IF  (TYPE. NE. PLATE)  GO  To  3000 


• - - . 


143 


r — 


c 

c 

c 

c 

c 

c 

c 

< 

c 

c 

c 


3UU0 


THEN  SET  BEAM-PLAT  £ FLAG  TO  INDICATE  PLATE  (2) 
3PFLAG  = 2 
GO  TO  9000 

ELSE  THIS  IS  A BAG  TYPEiTITLE  CAkO  — WRITE  MESSAGE 
wki.  T E (6,9901)  TYPE,  TITLE 

END  CASE  ON  TYPE 

END  IF  (EOF) 


9000  RETURN 


FORMAT  STATEMENTS 

9900  FORMAT  (16A5) 

9901  FORMAT  ( 1H1/1HG /1H0 ,*BAD  TYPE  & TITLE  CARD  FQUND:*/5X,  16A5/ 

1 * CARDS  MAY  BE  OUT  OF  ORDER*) 

ENO 


J 


14  4. 


bOi •- S. 


OOOO  OOOOO 


6UBRJJTINE  RHUAI A 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


THIS  SUBROUTINE  READS  DATA  CARDS  it  2 THRU  6 FOK  THE 
CJRKENT  CASE  AND  PRINTS  THE  INPUT  DATA  OUT  ON  THE 
LISTING.  IF  INPUT  ERRORS  OCCUR,  THEY  ARE  COUNTED  AND 
THE  NUMBER  OF  ERRORS  IS  PRINTED  AFTER  I HE  LISTING  OF 
THE  INPUT  VALUES.  IF  AN  END-OF-FIlE  IS  FOUND  , AN  ERROR 
MESSAGE  IS  PRINTED  AND  THE  EOF-OF-FILE  FLAG  IS  SET. 

THIS  ROUTINE  ALSO  USES  THE  NEWLY-KEAD  VALUES  OF  AIDA,  F, 

AND  WAVE -FUNCTION-CODE,  If  THEY  ARE  VALID,  T3  SET  UP 
SUBSCRIPTS  FUR  THE  TYPE  OF  CASE  DESCRIPTION  TO  QE  INCLUDED 
IN  THE  OUTPUT  PAGE  HEADINGS. 

ALL  INPUT  DATA  READ  IS  PLACED  IN  THE  COMMON  BLOCK,  INPUTS. 
THE  PRINT  HEADING  SUBSCRIPTS  ARE  PLACED  IN  THE  COMMON  BLOCK, 
PRINTS. 

THIS  SUBROUTINE  IS  CALLED  BY  THE  PROCEDURE , NXTCAS.  IT 
USES  THE  SYSTEM  SUEPROCE DURE , EOF,  TO  CHECK  FOR  AN 
ENU-OF-FILE  CONDITION.  IT  CALLS  THE  SU B PR 0 CE DURE , PAGE,  TO 
SET  UP  THE  HEADINGS  ON  A FRESH  PAGE  FOR  THE  NEW  CASE. 


INTEGER  DAJXFG, 
COMMON  /FLAGS  / 

1 

COMMON  /1NFUTS/ 

1 

2 

COMMON  /PRINTS/ 

1 


BPPLAG,  DON  EF  G , EOF  L AG 

bADXFG,  BPFLAG,  UO.NLFG,  EOFLAG,  I Ei\RF  G , KOUNT, 

LOADFG,  MECPFG,  MIST  FG,  NC  CNF  G 

A,  AIDA,  ALPHA,  AR,  BETA,  D,  F,  H,  PC,  PE, 

Q,  SIGMAC,  SIGMAR,  TAU,  TIN3R,  T M AX,  TPrInT  , 

XM,  WAVEFN 

FLAG,  MAXLIN,  NAiUA,  NF,  NUMLIN,  NUMPAG,  NWAVEF, 
STEFCT,  TIMNOW,  T I T L E ( 1 5) , TODAY,  TYPE 


SET  UP  OUTPUT  HE AO  I N GS  ON  A NEW  PAGE  AND 
INITIALIZE  COUNT  OF  INFUT  ERRURS  TG  LAST  COUNT 

NUMPAG  = C 
CALL  PAGE 
LCOUNT  = KOUNT 


READ  REMAINING  CARDS  FOR  CASE 
AND  PRI NT  OUT  INPUT  VALUES 


C 


READ  ( 5 , 991  0 ) 

1 

2 

3 

4 

white  (b,<mi> 

1 

2 

white  (6,9912) 
1 


A,  H,  AR,  XM, 

FC,  PE,  ALPHA,  TAU, 

F,  WAVEFN,  AIDA,  BETA, 
SIGMAC,  SIGMAR,  Q,  D, 

T I NCR,  TMAX,  TPRINT 
A,  H,  AR,  XM, 

PC,  PE,  ALPHA,  TAU, 

F,  WAVEFN  , AICA,  BE!  A 
SIGMAC,  SIGMAR,  G,  D, 
Tif.C.-v,  TMAX,  TPRINT 


14S. 


non  non  ooo  ooo  oooooo  ooo  o o o o n o ooo  no 


IF  end-of-file  FOUND 

IF  (.NOT.  (EOF 1 5 ) .NE.U.O)  ) GO  TO  1000 

THEN  SET  END-OF-FILE  FLAG  AND  PRINT  ERROR  MESSAGE 

EOFuAG  = 1 
WRITE  (6,9902) 

WRITE  (6,9904) 

ELSE  CONTINUE 
END  IF  (EOF) 

IF  ANY  INPUT  ERRORS  WERE  FOUND 

1000  NEKR  = KOUNT  - LCOUNT 

IF  (.NOT.  (NERK.NE.0)  ) GO  TO  2000 

THEN  SET  INPUT-ERROR  FLAG  AND  PRINT  ERROR  MESSAGE 

IERRFG  = 1 

WRITE  (6,9903)  NERR 

WRITE  (6,9904) 

ELSE  CONTINUE 
ENO  IF  (IERR) 

I-  AIuA,  F,  OR  WAVEFN  ARE  WITHIN  RANGE 

20  0 C IF  (.NOT.  (AIDA. GE. -1.0  .AND.  AIDA.  LE.  1.0)  .OR. 

1 .NOT.  (F.EQ.l.G  .OR.  F.EQ.2.0)  .OR. 

2 .NOT.  (WAVEFN. EQ.l. 0 .OR.  W A VEF N. EQ. 2 . 0 ) ) GO  TO  2500 

THEN  SET  UP  PRINTING  SU5SCRIPTS 

NAIOA  = AIQA  ♦ 2.0 
NF  - F 

NWAV6F  = WAVEFN 
GO  TO  4000 

ELSE  SET  ERROR  FLAG  AND  PRINT  ERROR  MESSAGE 

25  0 C IERRFG  = 1 

Wr.ITE  (6,9905) 

END  IF  (AIUA,  F,  WAVEFN) 

400C  KEIU  RN 

FORMAT  STATEMENTS 

9902  F U <v  M A T (1HL/ lHc,  *END  OF  FILE  FGUNs  DURING  I N ? J T OF  DATA*/ 

1 * CASE  IS  TERMINATED  AS  NOT  ALL  DATA  IS  PRESENT*) 

99 U 3 FORMAT  (1H0/1HO,  14,  * U AT  A EmROK(S)  wEkE  FOUND  IN  INPUT*/ 
1 * CASE  IS  TERMINATED  AS  DATA  Is  If)  INCORRECT  FuRM*) 


990  4 FORMAT 

1 

990  F FORMAT 
1 

9911  FORMAT 
9911  FORMAT 

1 6X 

2 7 X 


( 1HG, 

* 


nu , ♦SOME  OF  THE  VALUES  LISTEO  ABOVE  MAT  ACTUALLY  BE1 
- Fr\0 M PREVIOUS  CASE*) 

( 1HC/ 1H0 ,* VALUES  OF  Aj 
* CASE  IS  TERMINATED1 
' ' " ° n /)  , 3F12.0) 


OF  AIOA, 

* ) 


4 

5 
b 

7 

8 
9 
A 
B 
C 
0 
E 


CASE  __ 

( 4 ( 4 FI  2 . 0 
. (1HG, 

i 5 1 HPL  AT  E HALF  WIDTH  OR  BEAM 
.51HGEAM  OR  PLATE  THICKNESS 
, 5 1HLEN3  TH  TO  WIDTH  RATIO, 


F,  Ok  WAVEFN  ARE  OUT  OF  RANGE*  7 


• J * -A*  ' < U.  rA  I « V 

7X  > 5 1HLEN3  TH  TO  WIDT) 

7 X»  5 1H  M ASS  PER  UNIT  AREA,  LBS 
1 HO , 

bX ,51HMAXI MUM  DISTRIBUTED  PRESSURE 
7 X»51HUN IF  CRM  PRESSURE  LOAD,  PSI. 
7X,5lHPRESSURE  DECAY,  DIMENSIONLESS 
7Y.;^o3rccnyi:  hiid  attf.m.  SEC# 


HALF  SPAN,  IN 

IN. 

DIMENSIONLESS 


-SEC . SOD/IN. CUBED 
LOAD,  PSI. 


( A> , G15. 8/ 
TH)  , G15 . 8/ 
l AR) , 015. 8/ 
(XM) , G15. 8/ 

(PC) ,G15. 8/ 
(PE) ,Gl5. 8/ 
(ALPHA) , G15 . 8/ 

/ T a n -4  ir  a / 


(TAU) ,315. 8/ 


2 

3 

4 

5 
b 
7 

6 


7X,51HPRESSURE  DECAY,  DI 
7 X , 5 1H PRES  SUR E DURATION, 

1H0 , 

bX,51HSUP?0RT  FACTOR:  1=SIMPLY,  2=CLAMPEQ 

7X,51HWAVE  FUNCTION:  1=  GENERAL  , 2=SQUAR£  ,r,„.u.r„,, 

7X, 51HWEI3HT  VECTOR:  0 = VERT,  l = ExF  BLW , -1  = EXP  AB  V (AIuA), 

7X  ,51HSPAT  IAL  PRESSURE  DECAY  CONS T ANT , D I MENS  I CNlE SS ( 3E T A ) , 

« u A r / < un 

(SIGMAC) , G15 . 87 
( S I G M A R ) ,G15.8/ 
oc'  (0)  , G15. 87 


t iRL  r K L.  o JUK  t-  U tO  M T OUNS  I ANI  , Ui  lONLti 

9912  FORMAT  ( 1H0  , 

1 fcX , 51HCONCRE  T E COMPRESSIVE  STRENGTH,  PSI.  ( 

2 7X,51H RE IN  FORCED  STEEL  YIElD  STRESS,  PSI.  ( 

3 7X , 51HREIN FORCEMENT  RATIO  IN  TENSION,  DIMENSI ONLES 

4 7X» 5 lHREiN FORCING  DISTANCE,  IN. 

1H0, 


( F ) , G15 . 8/ 
(WAVEFN), G15. 8/ 
- 1 = E X P AB  V (AIuA)  ,G15.8/ 

i G 1 5 , 8 ) 


1H0, 

6X,51HTIME  INCREMENT,  SEC. 

7X.51HTIME  MAXIMUM,  SEC. 

7X,51HTIME  STEP  INTERVAL  FER  PRINTED  LINE, 


6X , 5 : 


SEC. 


(TINCR)  , G1 5 . 8 7 
(THAX) , G15. 8/ 
(TPRINT) , Gib. 8) 


(D) t G15 . 8/ 


End 
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SUDRuJTINE  JCNIRL  (NTRI.S) 


C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 


THIS  SUBROUTINE  STARTS  WITH  EA  - CASE  AT  TIME  T-j  AND 
AFTER  SOME  INITIALIZATION,  PUSHES  THE  CASE  THROUGH  THE 
TIME  STEPS,  USING  T INCH  AS  THE  TiMS  STEP  INCREMENT  AND 
TMAX,  AS  THE  fIME  STOP  INuICATuR. 

THE  BAD-HINGE-LOCATION-F LAG  IS  SET  IF  A 3AU  VALUE  FOR  THE 
HINGE  LOCATION  IS  DISCOVERED  BY  THE  SUBP RO C EG URE , CHEKXH * 
THE  lOGP-CONTRUL-FLAG  IS  SET  WHEN  THE  5UflP KOCE DURE , TSTEF, 
FINDS  THE  TIME-STEP-LOOP  SHOULD  BE  TERMINATED*  FUR  ANY 
REASON.  THE  CASE- IS-CUNE- FLAG  CAN  BE  SEl  IN  THE  FOLLOWING 
WAYS  : 

11  IF  THE  TIME-STEF-LOOF  HAS  oEE  N ATTEMPTED  MAXTRI 
TIMES  (SET  IN  TCNTRL) 

2l  IF  THE  ORIGINAL  CONSTANTS  FOR  THE  CASE  ARE  NEGATIVE 
OR  WOULD  NUT  CONVERGE  (SuT  IN  TZERD) 

3:  IF  THE  MAXIMUM  DEFLECTION  HAS  BEEN  FOUNC  OR  IF  THE 

TIME  MAXIMUM  HAS  bEEK  EXCEEDED  (SET  IN  TSTEP). 


THE  PARAMETER,  N TRIES,  IS  INCREMENTED  BY  THIS  ROUTINE, 
WHENEVER  THE  T 1ME-ST EF-LOOP  IS  TERMINATED,  TO  INDICATE  THE 
NUMBER  OF  TIMES  THE  LOOP  HAS  BEEN  TkIEO.  IN  THE  CASE  OF 
AN  UNSUCCESSFUL  TERMINATION,  THE  LOOP  WILL  BE  TRIED  AGAIN, 
UP  TO  MAXTRI  TIMES,  WITH  THE  TIME  STEF  HALVED  EACH  TIME. 


THIS  SUBROUTINE  CALLS  THE  FOLLOWING  SUSP RO CEOUKES i 

TZ ERG!  TO  INITIALIZE  CONSTANTS  AND  VARIABLES  FOR  THE  CASE 
TSTEF:  TO  COMPUTE  THE  DESIRED  VARIABLES  AT  EACH  TIME  STEP 

FROM  1'IME  = T I NCR  THROUGH  THE  TIME,  T = TMAX. 


THIS  RCUTINE  IS  CALlED  BY  THE  PROCEDURE,  DRIVER. 


integer  babxfg, 

COMMON  /FLAGS  / 

1 

COMMON  /INPUTS/ 

1 

? 

COMMON  /RESULT/ 

1 

COMMON  /PRINTS/ 

1 


BPFLAG,  CCNEFG , EuFLAG 

BADXFG,  BPFlAG,  DONEFG,  EOFLAG,  IERRFG,  KOUNT, 

LOA  CFG,  MECPFG,  M1STFG,  NCO.NFG 

A,  AIDA,  ALPHA,  AR,  BETA,  D,  F,  H,  PC,  PE, 

Q,  S1GMAC,  SIGMAR,  TAU,  T1 NCR , T MAX,  T PRINT , 

XM,  WAVEFN 

DELTA,  T,  THETA,  THETAD,  V,  VtL,  WF,  WK,  WP , 

X,  XH 

FLAG,  MAXLIN,  NAIDA,  NF,  NUMLlN,  NUMPAG,  NWAVEF, 
STEPCT,  TIMNDK , T1TLE(15>,  TODAY,  TYPE 


DATA  MAXTRI  /2/ 
DATA  BLANK  /1H  / 


SET  MECHANISM-FLAG  TO  USE  MECHANISM  2 TO  CHECK  XHL , 

clea-<  bad-hinge-location  and  loof-ccntfol  flags, 

SET  THE  FIRST -TIME-FLAG  FOR  THE  FIRST  TWO  PASSES  THkU  CUMP, 
CLEAR  NON-CONVEKGEnT-FL AG  FOR  PLATE  CASE, 

CLEAR  FAILURE  FLAG  FOR  PRINTOUT, 

AND  SET  TYPE-GF-LCAD  FLAo  TU  BLAST  LOAD  UNLESS  UETA=0. 

also,  initialize  time  variables,  compute  case  constants 

146. 


. 


....... 


: FROM  INPUT  AND  PRINT  T = J kESULTb  ON  NEW  PAGE  ( T ZERO ) • 

MECHFG  = 2 
3 A ()  X F 3 = U 
L OC  PF 3 = 0 
H1STFG  = 2 
NCONFG  = 0 
L0A0F3  = 2 

IF  (BETA. EG.  0.  G ) LOADFG  = 1 

r 

THETA 
TmETAJ 
BELT  A 
7 EL 
WF 
Wp 
\J 

WK 

FLAG  = BLANK 
NUMLIN  = MAXLIN 
CALL  TZEkO  (NTRIES) 

IF  THE  ORIGINAL  CONSTANTS  WERE  CALCULATED  SUCCESSFULLY 
IF  (OONEFG.EQ.l)  GO  TO  5000 
THEN 

DO  WHILE  lucp-comrcl-flag  IS  clear 

(I.E.,  WHILE  TIME<TIHE  MAX  AND  BADXFG  IS  CLtAR) 

1000  IF  (LUCPFG.NE.G)  GO  TO  4000 

CALL  TSTEP  <LOOPFG) 

GO  TO  1000 

END  WHILE 

ADO  1 TO  NUMBER  OF  TRIES 
4 0 0 C NTRIES  = NTRIES  «■  1 

IF  This  WAS  THE  LAST  TRY  ALLOWED,  SET  CASE-IS-D ON E-FLAG 
IF  (NTRIES. GE.MAXTRI)  LONEFG  = 1 

ELSE  CONTINUE,  SKIPPING  CALCULATIONS 
END  IF  (ORIGINAL  CONSTANTS) 


= 0.0 
= 0.0 
= 0.J 
= 0.0 
= 0.0 
= 0.0 
= 0.0 
= 0.0 
= 0.0 


5 0 0 j Kt TURN 


o o o o o o o 


SU8RGJT  INE  T/EPO  (NTRIES) 


THIS  SU3R0UT  I lit  COMPUTES  AMU  PRINTS  CA^E  CONSTANTS  PROM  THE 
GIVEN  INPUT  VALUES.  IT  ALEC  PROVIDES  FOR  THE  FIRsT  lINE  OF 
OJTPUT  ( T = 0 ) TO  DE  PkINTEC  OUT  ON  A NEW  PAGE  WITH  HEADINGS. 

'THIS  SUBROUTINE  IS  CALLED  BY  THE  PROCEDURE  TCNTRL  AND  CALLS 
THE  FOLLOWING  SUBF ROCEDUKES : 

3T  ZERO l TO  INITIALIZE  BEAM  CONSTANTS  AND  VARIABLES 
PT  ZERO  : TO  INITIALIZE  PLATE  CONSTANTS  ANO  VARIABLES 
CALXHO:  TO  FIND  THE  OkIC-INAl  HINGE  LOCATION,  XHQ 

CHEKXH!  TO  CHECK  THE  OF.IGINAL  LOCATION  OF  THE  HINGE , 

XHG,  AT  TIME  = G 

PRINTRI  TO  PRINT  OUT  THE  INITIAL  VALUES  FOk  THE  GIVEN 
CASE  AT  TIME  = 0 

EXP!  THE  SYSTEM  LIBRARY  FUNCTIUN  FOR  EXPONENTIALS 

THE  INPUT  PARAMETER,  NTRIES,  IS  USED  TO  PREVENT  THE  CALCU- 
LATION ANO  PRINTING  OP  THE  CCMPUTEC  CONSTANTS  FOR  THE  CASE 
ON  THE  SECOND  AND  SUCCEEDING  TRIES.  THE  PARAMETER  WILL 
NOT  BE  ALTERED  IN  ANY  WAY. 

THE  CASE-I S-OONE-FLAG  MAY  BE  SET  IN  TWD  WAYS!  FIRST, 

WHEN  THE  PLATE  CONSTANT,  Z,  CANNOT  BE  FOUND  OR  GOES  NOT 
CONVERGE,  AND  SECOND,  IF  NEGATIVE  CONSTANTS  ARE  COMPUTED. 


integer  badxfg,  bpflag,  bcnefg,  eoflag 

COMMON  /FLAGS  / BADXFG,  BPFLAG,  DCNEFG,  EOFLAG,  1ERRFG,  KOUNT, 
1 LOADFG,  MECHFG,  M1STFG,  NCONFG 

COMMON  /INPUTS/  A,  AIDA,  ALPHA,  AR,  BETA,  D,  P,  H,  PC,  PE, 

1 Q,  SIGMAC,  SIGMAk,  TAU,  TINCR,  TMAX,  TPRINT, 

Z XM,  WAVEFN 

COMMON  /RESULT/  DELTA,  T,  THETA,  THETAD,  V,  VEL,  WF,  WK,  WP, 

1 X , XH 

COMMON  /CONSTS/  ACUjE,  ARSU,  ARSZP1,  ARZSP1,  ASQ,  B,  BARA1, 

1 BAkCI,  QAkGI,  BARH1,  6ARJ1,  BETACB,  BE  T AFP , 

2 BET ASQ,  DELTAK,  EPSLNU,  EXP3ET , FOURTH,  HALF, 

3 ONEMZ , GNEPZ,  SIXTH,  fHETUl,  THIRD,  W,  XMU,  Z, 

4 ZCUBE,  ZFOUR » ZSQ 


IF  THIS  IS  THE  FIRST  TIME  THkU  TZEKJ  FOR  THIS  CASE 

IF  (NTRI ES.NE. C)  GO  TO  3000 

THEN  USE  INPUT  VALUES  TO  COMPUTE  CONSTANTS 

XMU  = V.U  * 0**2  * Q * SIGMAk  * 

1 (1.0  - 0.59  * Q * (S1GMAR/SIGMAC)  ) 

AiQ  = A * A 
ACU3E  = A**3 
ARuO  = AF*AR 
BETASQ  = JET  A*HET  A 
BETACB  = DETA**3 

ISO  . 


GETAFR  = 3 ET  A **  A 
EX POET  = EXP(-DETA) 

W = 386. A * XM 

OE  L f A<  = AIDA  * 386. A 

T HE  r U 1 = ( A. C * H * EPSLNU)  / (u*"2) 

C 

C CASE  ON  TYPE  OF  CASE 

C 

C IF  BEAM  CASE 

C 

IF  (3PFLAG.EQ.1)  CALL  BTZERC 
C 

C IF  PLATE  CASE 

C 

IF  ( oPFLAG.EQ. 2)  CALL  FTZERG 
C 

C END  CASE  ON  TYPE 

C 

C IF  Z CALCJLATIONS  DID  NOT  CONVERGE 

C 

IF  (NCONFG.NE.l)  GO  TO  1000 
C 

C THEN  SET  CASE-IS-CONE-FLAG  TO  TERMINATE  CASE 

C 

DONEFG  = 1 
GO  TO  3000 
C 

C ELSE  FIND  VALUE  OF  ORIGINAL  HINGE  LOCATION 

C AND  PRINT  OUT  COMPUTED  CONSTANTS 

C 

1000  CALL  CALXH0  (X0,  XH0) 

C 

IF  (BPFLAG. EQ. 1)  WRITE  (6,9901)  B,  XMU,  W , XH0 
IF  (BPFLAG. EQ. 2)  WRITE  (6,9902)  3,  Z,  XMU,  W,  XH0 
C 

C IF  ANY  NEGATIVE  CONSTANTS  WERE  COMPUTED 

C 

IF  (.NOT.  (B  .LE.0.0  .OR. 

1 XMU. LE.0.0  .OR. 

2 W .LE.0. 0 ) ) GO  TO  3000 

C 

C THEN  PRINT  ERROR  MESSAGE  AND  TERMINATE  CASE 

C 

WRITE  (6,9903) 

UONEFG  = 1 
C 

C ELSE  CONTINUE 

C END  IF  (NEGATIVE  CONSTANTS) 

C 

C ENU  IF  (Z  CONVERGED) 

C 

C ELSE  CONTINUE 

C ENU  IF  (1ST  TIME) 

C 

C Ir  CASE  CONSTANTS  APPEAk  OKAY 


151. 


Oooo  oooo  ooo 


3  0 0 G IF  (JONEF&.NE. 0)  GO  TO  4000 

THEM  INITIALIZE  X AND  XH  TO  ORIGINAL  HINGE  LOCATION* 

CHECK  MECHANISM  AND  PRINT  TIME  ZERO  RESULTS  CN  NEW  PAGE 

X = X0 
XH  = XHG 
CALw  CHEKXH 
CAL..  PRINT  R 

ELSE  CONTINUE 
END  IE  (CONST  A N'T  OKAY) 

40 00  RETURN 

FORMAT  STATEMENTS 

9901  FORMAT  ( 1 HQ/  1HQ,  2X,  ^COMPUTED  CONSTANT  VALUES*  / 1H0, 

1 6X  ,5 1HBEAM  WIDTH,  IN.  (61,015.8/ 

2 ZX , 5 1HHINGE  MOMENT,  IN. -LBS. /IN.  ( X MU > , G15 . 8/ 

3 ZX, 51HWEIGHT  PER  UNIT  AREA,  L6S./IN.SQ.  (W),G15.8/ 

4 ZX,51H0RIGINAL  HINGE  LOCATION,  IN.  <XH),G15.8) 

9902  FORMAT  ( 1HC/  1H0,  2X,  *COMPUTEC  CONSTANT  VALUES*/  1H0, 

1 6X.51HPLATE  HALF  LENGTH,  IN.  <3),G15.8/ 

2 ZX,51H RATIO  OF  FINAL  HINGE  LCC  TO  &,  Ol ME NSI ONl E S S (Z), 015.8/ 

3 ZX.51HH1NGE  MOMENT,  IN. -LBS. /IN.  (XMU),Gl5.8/ 

4,  ZX,51HWE1GHT  PER  UNIT  AREA,  LBS. /IN. S3.  (K),&13.8/ 

5 ZX, 51HORIG INAL  HINGE  LOCATION,  IN.  (XH),&15.8) 

9902  FORMAT  ( 1H0/  1H0,  *SOME  COMPUTED  CONSTANTS  FOR  THIS  CASE  * 

1 * ARE  NEGATIVE.*/  IX,  *SGME  INPUT  VALUES  MUST  BE  IN  ERROR.*/ 

2 IX,  * C ASE  IS  TERMINATED.*/) 

C 

End 


15  2 


— 1 ***— ~—»' w 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 


c 

c 

c 

c 

c 

c 


c 

c 

c 
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SUBROUTINE  arZEkO 

THIS  SUBROUTINE  CALCULATES  THE  INITIAL  VALUES  OF  CONSTANTS 
AND  VARIABLES  TO  BE  USED  BY  THE  GIVEN  BEAM  CASE. 

ALL  OF  THE  CONSTANTS  AND  VARIABLES  INITIALIZED  ARE  PASSED 
TO  THE  OTHER  PROCEDURES  THROUGH  THE  COMMON  BLOCKS,  COMPS 
AND  CONSTS. 

THIS  SUBROUTINE  IS  CALLED  BY  THE  PROCEDURE,  T ZERO.  IT  CALL 
NO  SUB  PROCEDURES. 


INTEGE 

R BAOXFG, 

BPFLAG,  DON  EFG , E OF L AG 

COMMON 

/flags  / 

OADXFG,  GFFLAG,  DjNEFG,  EOFLAG,  1ERRFG,  KCUNT, 
LOADFG,  MECHFG,  M1STFG,  NCON^G 

COMMON 

/COMPS  / 

ATHE3B,  ATHEDL,  ATHEC1,  ATHED2,  AVDOT,  AWFub, 
AWFOL,  BTHED2,  31HED3,  BVOOT,  3WFDOT,  BWPDOT, 
PWFDB , PWFDL,  FWFDOT,  PVDOT 

COMMON 

/INPUTS/ 

A,  AIDA,  ALPHA,  AR,  BETA,  D,  F,  H,  PC,  PE, 

Q,  SIGMAC,  SIGMAR,  TAU,  T1NCR,  IMAX,  TPRINT , 
XM,  WAVEFN 

COMMON 

/CONSTS/ 

ACUBE,  ARSQ,  ARSZF1,  ARZSP1,  ASQ,  3,  BARA1, 
BARC1 , BARG1,  BARH1,  BARJ1 , 3ETACB,  BETAFk, 

BET  AS  G,  DELTAK,  EPSLNU,  EXP3ET,  FOURTH,  HALF, 
ONEMZ , CNEPZ , SIXTH,  THE TUI , THIRD,  W,  XMU , Z, 
ZCUBE,  ZFOUR , ZSQ 

SET  UP  BEAM  CONSTANTS 


3 = 2.  0 * A 

/ 

AR 

8 THE  02  = 3.0 

¥ 

HALF  * 

DELTAK 

/ 

A 

BTHED3  = 3.0 

¥ 

F * XMJ 

/ (XM 

¥ 

ACUBE) 

BWFOOI  = 2.0 

¥ 

ASQ  * B 

BWFDQT  = 2.0 

¥ 

G * F * 

XMU 

BVDOT  = 2.0 

¥ 

ASQ  * B 

* AIDA 

¥ 

W 

AT  hEu2  = BTHED2  8THED3 
AVDOT  = BVDOT  * HALF 

IF  THIS  IS  A LINEAR  LOAD 

IF  (LOADFG.NE. 1 ) GO  TO  20G0 

THEN  DETERMINE  LINEAR  LOAD  CuNSTANTS  FOR  COMP 

AT  HE  DL  = (3.u*PE  ♦ 2.0*PC)  / (2.0  * A * XM> 

AWFOL  = BWFDDT  * (PE*hALF  + FC*‘THxRO) 

GO  TO  30  C 0 

ELSE  DETERMINE  BLAST  LOAD  CONSTANTS  FOR  COMP 

ATHEDJ  = < 3.0/ (XM+A)  >MPE*HALF 

♦ (PC/dET ACo)  MbtTASQ  - 2.0*C£TA  f 2.0 

- 2. u*  EXPBET ) ) 


1S3  . 


oo  non  oooo 


SUBROUTINE  PTZERO 


:: 


* 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


THIS  SUBROUTINE  COMPUTES  INITIAL  VALUES  OP  VARIABLES  AN'O 
CONSTANTS  TO  BE  USED  bY  THE  GIVEN  PLATE  CASE.  THE  PLATE 
CONSTANT  Z IS  FOUNL  USING  THE  ROOT-FINDING  ROUTINE,  BISECT, 
WITH  Th£  Z FUNCTION,  FTNZ . A MESSAGE  IS  PRINTED  AND  THE 
CASE  IS  TERMINATED  IF  Z CANNOT  bE  FOUND. 

THIS  SUBROUTINE  IS  CALLEO  BY  THE  PROCEDURE,  T ZERO,  AND 
CALLS  THE  FOLLOWING  SUBPROC E DURE S I 

BISECT:  TO  SOLVE  FCk  Z (THE  FUNCTION,  FTNZ , IS  CARRIED  TO 

BISECT  AS  A PARAMETER) 

CALOARI  TO  SET  UP  EQUATION  OF  MOTION  CONSTANTS  FOR 
MECHANISM  1. 

THE  NGN-CONVERGENT-FLAG  IS  RETURNED  TO  THIS  ROUTINE  BY 
THE  ROUTINE,  BISECT,  AS  SET,  IF  THE  ROOT  FOR  THE  Z FUNCTION 
CANNOT  BE  FOUND.  THIS  IS  PASSED  THROUGH  THE  COMMON  BLOCK, 

Fi_  AGS , BACK  IG  THE  CALLING  PRuCELURE,  T ZERO  • 

ALL  THE  CONSTANTS  ANC  VARIABLES  INITIALIZED  ARE  STORED  AND 
PASSED  TO  THE  OTHER  PROCEDURES  BY  THE  COMMON  BlOCKS,  BARS, 
COMPS,  AND  CONSTS. 


E X T E RN  A L F T NZ 

INTEGER  8A0XFG,  BP  FLA  G,  CONEFG,  EOFLAG 


COMMON 

/FLAGS 

/ 

EADXFG,  BPFLAG,  DONEFG,  EOFLAG 
LOADFG,  MECHFG,  M1STFG,  NCONFG 

COMMON 

/BARS 

/ 

BARA,  BARB,  OARC,  BARD,  BARDEL 
BAKH,  BARB,  BARK,  BARK,  BARS, 
EXP3X1,  EXP  B ZX 

COMMON 

/COMPS 

/ 

ATHEDB,  ATHEDL,  ATHED1,  ATHE02 
A WFOL , BTHE0  2,  BTHED3.  BVDOT, 
PWFOB,  F WFD  u , PKFDOT,  PVDOT 

COMMON 

/INPUTS/ 

A,  AIDA,  ALPHA,  AR,  BETA,  D,  F 
U , SIGMAC,  SIGMAR,  TAU,  TINCR, 
XM,  WAV  EFN 

, IERRFG,  KOUNT, 

, BARCNM,  BARG, 
BART,  BARTHE,  BA \U 

, AVQGT,  AWFCB, 
8WFD0T,  6WPDOT, 

, H,  PC,  PE, 

TMAX,  TPRINT , 


COMMON  /CONSTS/ 

1 

2 

3 


ACUuE,  ARSQ,  ARSZf-1,  ARZSP1,  A SQ,  B,  BARAi, 
BARC1,  BARG  1 , 3ARH1,  3ARU1,  3ETAC6,  8ETAFR, 

BET AS  Q,  DELTAK,  EFSLNU,  EXPBET,  FOURTH,  HALF, 
ONENZ,  GNEPZ , SIXTH,  THETU1,  THIRD,  W,  XMU,  Z, 
ZCUt3E  , ZFGUR,  ZGQ 


31 


i 


SOLVE  FOR  Z 


NCCNF3  = 0 

CALL  3ISEC1  (-0.0001,  l.JOCi,  Z,  NCONFG,  FTNZ) 
IF  Z IS  NUI  IN  THE  INTERVAL  (0,1) 

IF  (NCONFG. EQ. 0 ) GO  TO  1200 

THEN  PRINT  MESSAGE  AND  TFRM1NATE  CASE 

155. 


4 


RkITE  (6,9901) 
GO  TO  3GC0 


ELSE  CONTINUE  UN 


e 

A * AR 

ZSQ 

Z * Z 

ZCU3E  = 

Z**3 

ZFOJR  = 

Z**  4 

ARSZF1  = 

ARSQ 

* Z ♦ 

1.0 

ARZSPl  = 

AR  * 

ZSQ  + 

1 . 0 

ONEMZ  = 

1.0- 

z 

ONEP  Z = 

1.0  + 

z 

BARA  1 = 

1.  0 / 

(ZSQ  * 

BEdF  A 

* XM  * 

A) 

B ARC  1 = 

- DELTA*  / A 

BARG1  = 

- <F  * 

XMU ) 

/ (XM  * 

ACUBE  * 

AR) 

B ARM  1 = 

1.  0 ♦ 

AR 

BARJ1  = 

1.0/ 

( XM  * 

BET  ACS) 

CAL^  CALBARd. 

0) 

PWPDOT  = 

4.0  * 

F * XMU  * A * 

ARSZP1 

/ (Z  ♦ A 

PVOOT  = 

4.0  * 

A1UA  * 

W ♦ ASQ 

* B 

ATHEDl  = 

BART  HE 

AT  HE  02  = 

B ARC 

* BARD 

+ BARG  * BARH 

AVDQT  = 

P 0 DOT 

* (1.0 

- ONE PI 

♦HALF  + 

Z*  T HI RQ) 

IF  THIS  IS  A LINEAk  LOAD 

IF  (LOADFG .NE.l)  GO  TO  2000 

THEN  DETERMINE  LINEAR  LOAD  CONSTANTS  FOk  COMP 

PRFDL  = 4.0  * AR  * ACUBE 
AHFOL  = FWFOL 

* (PE  ♦ (1.0  - CNEPZ¥HALF  + Z* THIRD) 

♦ PC  * (THIRD  - CNEPZ*SIXTH  + SIXTH  - Z * ONEP Z * FOURTH ) ) 
GO  TO  3000 

ELSE  DETERMINE  BLAST  LOAD  CONSTANTS  FOR  COMP 

PHF03  = 4.0  * ASQ  * B / (ZSU  * BETAFR) 

ANFOb  = FWFDd 

* (PE*bETAFR*ZSQ*  (1.0  - Z’HAlF  + ONEPZMHIRD) 

♦ PCM2.0*BETA3Q*ZSQ  - 4.0*BETA*ZSQ 

♦ EXPBZX*  ( - 3.0*BETA*Z  ♦ Z.O^BETA+Z  ♦ 3.0) 

♦ EXP3X1*ZSQ* (BET ASQ  - 5.0’OETA 

- BET  ASQ  *■  2.0*BETA  ♦ 9.0) 

- EXPbET*  ( 2. 0*BETA*Z*GNEPZ  ♦ 9.0*ZSQ  3.0)  ) ) 

END  IF  (LINEAR  LOAD) 


END  IF  (BAD  Z> 


SUBROUTINE  CALXHO  (XG,  XHO) 


THIS  SUBROUTINE  HANDLES  CALLING  THE  METHOD  FOR  SOLVING 

whatever  function  is  necessary  to  find  the  original  hinge 

LOCATION  FOR  THE  CASE--OR  PRINTING  AN  ERROR  MESSAGE  IF 
IT  CANNOT  BE  FOUND  EXACTLY. 


THIS  ROJTINE  IS  CALLED  BY  1HE  SUBROUTINE,  TZEkO. 

THE  ONLY  SUBPROCEDURE  CALLED  BY  THIS  ROUTINE  IS  THE 
ROOT-FINDING  SUBPROCEDURE,  BISECT,  WHICH  USES  THE 
FUNCTION,  BFTNX  OR  PFT NX , DEFENDING  OF  THE  TYPE  OF  CASE. 

THE  PARAMETERS  RELATING  TO  THE  COMPUTED  HINGE  LOCATION, 

XO  AND  XHO,  ARE  RETURNED  TO  THE  CALLING  PROCEDURE,  T ZERO . 


EXTERNAL  BFTNX,  PFT NX 

INTEGER  BAOXFG,  BPFLAG,  LUNEFG,  EOFLAG 

COMMON  /FLAGS  / BADXFG,  BPFLAG,  DONEFG,  EOFLAG,  IERKFG,  KOUNT, 
L LOADFG,  MECHFG,  M1STFG,  NCONFG 

COMMON  /INFUTS/  A,  AIDA,  ALPHA,  AR,  BETA,  D,  F,  H,  PC,  PE, 

Q,  SIGMAC,  SIGMAR,  TAU,  TINCR,  f M AX , TPRINT , 

! XM,  WAVEFN 


: COMPUTE  XHO  IN  THE  INTERVAL  ( C, A) 

■\ 

» 

DETERMINE  TYPE  OF  CASE  FOR  CORRECT  XHO  FUNCTION 

/ 

) IF  BEAM  CASE 

IF  (BPFLAG. EQ.  1) 

1 CALL  BISECT  (0.00001,  0. 90999,  XO,  NCONFG,  BFTNX) 

/ 

; i - plate  case 

IF  (BPFLAG.  EQ.2) 

1 CALL  BISECT  (0.00001,  0.99999,  XG,  NCONFG,  PFTNX ) 

/ 

i END  CASE  OF  TYFE  OF  CASE  FOR  XHO  CALCULATION 

XHO  = XO  * A 

; CHECK  CONVERGENCE 

; CASE  II  BISECT  DID  NOT  FIND  A ROOT  IN  THE  INTERVAL  ( 0 , A ) 

/ 

IF  (NCONFG. NE. 1)  GO  TO  2100 
; ASSUME  XHO  = A 

WRITE  (6>»  990  3) 

NCONFG  = 0 
XG  = 1.0 
XHO  = A 

CASE  21  BISECT  DID  NOT  CONVERGE  AFTER  100  ITERATIONS 
2 1 0 L IF  (NCONFG. NE.-l)  GO  TO  3000 

158  . 


> * 


ASSUME  ROOT  FOUND  IS  CORRECT 
WRITE  (6,9904)  XHO 
NCONFG  = 0 

END  CASE  ON  CONVERGENCE 

3000  RETURN 

FORMAT  STATEMENTS 

9903  FORMAT  (/1H0,*NQ  VALUE  OF  THE  ORIGINAL  HINGE  LOCATION,  XHO,* 

1 * WAS  FOUND  IN  THE  INTERVAL  ( 0, A ) . */ 

2 * IT  IS  ASSUMED  TO  BE  THE  VALUE  OF  A.*) 

9904  FORMAT  (/1HQ,*THE  BISECTION  METHOD  USED  TO  FIND  THE  ORIGINAL* 

1 * HINGE  LOCATION,  XHO,  DID  NOT  CONVERGE  AFTER  100* 

2 * ITERATIONS.*/ 

3 * THE  RESULT  OF  THE  LAST  ITERATION  WILL  BE  ASSUMED* 

4 * CORRECT  « XHO  = *,  G15.8,  *.*) 


n o o o o 


SUBROUTINE  B IS  E Cl  (GUESSl,  GUESS2,  <001,  ERkrLG,  FTN) 


THIS  SUBROUTINE  TRIES  TO  FIND  A ZERO  ROOT  FO*  THE  GIVEN 
FJNCTION  WITHIN  THE  INTERVAL  (oUESSl,  GUESS  2)  USING  THE 
BISECTION  METHOD.  GUESS  1 MUST  BE  LESS  THAN  GUESS  2. 
NOTE:  THIS  kGUTINE  MAY  BLOW  uP  IF  THE  FJNCTICN,  FTM,  HAS 

A DISCONTINUITY  IN  THE  GIVEN  INTERVAL. 


THE  INPUT  PARAMETERS  ARE  DEFINED  AS  FOLLOWS: 

GUESS  1 » THE  LOWER  END  OF  ThE  INTERVAL  BEING  CONSIDERED 
GUESS2J  THE  UPPER  ENJ  CF  THE  INTERVAL  BEING  CONSIDERED 
FTN:  THE  EXTERNAL  FUNC11UN  SUBPROCEDURE  WHICH  COMPUTES 

THE  FUNCTION  FCR  WHICH  BISECT  IS  ATTEMPTING  TO 
TO  FIND  A ROOT. 


THE  OUTPUT  PARAMETERS  ARE  DEFINED  AS  FOLLOWS: 


ROOT: 


ERRFLGi 


THE  ROCT  FCR  THE  FUNCTION,  FIN,  IN  THE  GIVEN 
INTERVAL  , IF  FOUND 

THE  MIDPOINT  OF  THE  INTERVAL,  IF  THE  ROOT  WAS 
NOT  FOUND 

0,  THE  ROOT  WAS  FOUNO  SUCCESSFULLY 

-1,  THE  ROOT  FAILED  TO  CONVERGE  AFTER  100 

ITERAT IONS 

+ 1,  THE  FOOT  COULD  NCT  BE  F J UNO  AS  FTN(GUESSl) 
AND  FT  N ( GUESS2 ) HAVE  THE  SAME  SIGN,  INFERRING 
THAT  THERE  IS  EITHER  NO  ROOT  OK  MORE  THAN  ONE 
ROOT  FOR  THE  FUNCTION  IN  THIS  INTERVAL. 


THIS  ROUTINE  IS  CALLED  BY  THE  SUBPROCEDURES , FTZEkO  AMD 
CALXHO.  IT  CALLS  THE  SU3PR0CELURE,  FTN , WHICH  MAY  BE  ONE 
OF  THE  FOLLOWING  FUNCTIONS: 

FTNZI  CALLEU  FROM  PTZERO,  TO  FIND  Z 

OFTNX:  CALLED  FROM  CALXHO,  TG  FIND  XHO  FOR  A BEAM  CASE 

PFTNX  t CALLED  FROM  CALXHO,  TO  FIND  XHO  FOR  A PLATE  CASE. 


INTEGER  ERRFLG 
EXTERNAL  FTN 


CLEAR  ERROR-FLAG  AND  INITIALIZE  END  POINT  VALUES 
INITIALIZE  EPSILON  AND  PUT  UUMMY  VALUE  IN  F3 


ERRFLG  = 0 

ITER  = 0 

XI  = S UE  SSI 

X 2 — *jUlSS2 

EPSILN  = AbS ( X 2 - XI 

F 5 = 1.0 

SF  3 = SIGN  (1.0,  F 3 ) 


XI)  *•  0.C0JJ  0 0 01 


Fl  = FTN  C A 1 ) 

F2  = FTN  <X2  ) 

SF  1 = SIoN  ( 1.  0 , F 1) 
SF2  = SIbN  ( 1.  0 , F 2 ) 
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l 


I 

c 

C I-  THE  SIGNS  OF  Fl  A NL  F2  ARE  THE  SAME 

I C 

IF  ISF1.NE.SF2)  GU  TO  1000 

THEN  THERE  IS  AN  EVEN  NUMBER  OF  ROOTS  IN  THE  GIVEN  INTERVA. 
(0,2, 4, . . .) --SET  THE  ERROR  FLAG 

ERRFLG  = 1 

ROOT  = (X2  - XD/2.0  + XI 
GO  TO  5000 

ELSE  CONTINUE 

JO  WHIlE  F3  IS  NOT  ZERO  AND  ITER  < 100 
AND  (X2  - XI)  IS  NOT  ZERO 

1 0 0 G IF  ( ( SF3*F3)  .LT.  EPSILN  .OR. 

1 ( X 2 - XI)  .LT.  EPSILN  ) GO  TO  490  0 

ITER  = ITER  +■  1 

IF  (ITER  .GT.  100)  GO  TO  4800 

FIND  X3  BETWEEN  XI  AND  X2.  COMPUTE  F3. 


X3  = XI  + ( X2  - XlJ/2.0 
F3  = FTN  C X 3 ) 

Sr  3 = SIGN  (1.0,  F 3) 

IF  NO  kCGT  BETWEEN  XI  AND  X3 

IF  (SF3.NE.SF1)  GO  TO  2000 

THEN  MOVE  XI  TO  X3 

XI  = X 3 
Fl  = F 3 
SFl  = SF  3 
GO  TO  3000 

else  move  x2  to  X3 

20  00  X 2 = X 3 

F2  = F 3 
SF2  = SF3 

ENO  IF 

3000  GJ  TO  1000 

ENO  WHILE 

TAKE  CARE  OF  ERkOk — DIG  NCT  CONVERGE  IN  100  TRIES 
48  0 C ERKFLG  = -1 
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FUNCTION  FTN  7 (TRIALZ) 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 


THIS  FUNCTION  CALCULATES  A VAt_UE  FOR  T HE.  PL  Aft  Z FUNCTION 
FOR  THE  GIVEN  TRIAL  Z.  IF  TRIAL  Z IS  A CORRECT  GUESS 
AT  Z,  THE  kESULT  WILL  BE  IE KC. 

THIS  FUNCTION  IS  CALLED  BY  THE  ROOT-FINDING  S UBPROCEOURE , 
BISECT,  WHICH  WAS  CALLED  BY  THE  SU3PR0CEDURE,  F T ZERO • IT 
CALLS  NO  SUBPROCEOURES  OTHER  THAN  THE  SYSTEM  LIBRARY 
ROUTINE,  EXP. 


1 

1 

2 

3 

4 

1 

2 


INTEGER  tiADXFG, 
COMMON  /FLAGS  / 

COMMON  /CONSTS/ 


COMMON  /INPUTS/ 


BPFLAG,  DONEFG,  EOFLAG 

BADXFG,  BPFLAG,  DONEFG,  EOFLAG,  IERRFG,  KOUNT 
LOAOFG,  MECHFG,  M1STFG,  NCONFG 
ACUOE,  ARSQ,  ARSZP'l,  ARZSPl,  A SO , B,  BARA1, 
BARCl,  BARG1 , BARH1,  3 AR Jl , 3ETACB,  BETAFR, 
BET  AS  Q,  DELTAK,  EFSLNU,  EXP  BET  , FOURTH,  HALF, 
ONEMZ , ONEPZ,  SIXTH,  THETU1 , THIRD,  W,  XMU,  Z 
Z CUBE , ZFCUR,  ZSQ 

A,  AIDA,  ALPHA,  AR,  BETA,  D,  F,  H,  PC,  PE, 

Q,  SIGMAC,  SIGMAR,  TAU,  T I NCR , T M A X , TPRINT  , 
XM,  WAVEFN 


TRZSQ  = TRIALZ  * TR1ALZ 

IF  THIS  IS  A LINEAR  LOAO 


IF  (lOAUFG.NE.1)  GO  TO  1000 


THEN  USE  LINEAR  LOAD  FORMULA  FOR  Z 


FT NZ  = PE* (4.0*ARSG* TRIALZ  * 8.0*7KlA_Z  -12.0) 

1 ♦ PC* (2.0*  ARSQ*  TR1ALZ**3  + 3.0*TRZSQ  - 5.0) 

GO  TO  5000 

ELSE  USE  BLAST  LOAD  FORMULA  F Or.  Z 

IlOL  EXP3Z1  = EXP  (BETA  * (TRIALZ  - 1.0)  ) 

FT  NZ  = AF.SCl*TRIALZ*(FE*BETAFf  *TRZSQ*TRlALZ*SIXrH 
1 + PC*(EXF3Z1*(6..J  - 4 ,0*8ET  A*TRIALZ 

Z *BET  ASQ*  T RZ  SO ) 

3 - EXFBETMb.O  + 2 . 0 * 3 ET A * TRI AL Z ) ) ) 

4 ♦ PE*BET A FR*TRZSQ*( TRIALZ* THIRD  - hAlF) 

5 * PC*  (EXPBZ1*  (3.0  - 3. 0 *6ET  A*TrUAcZ  ♦ at  T ASQ*TRZSQ> 

5 * EXPBET* (2.0*3 ETA*TRZSO  ♦ 9.0*TRZSC  - 3.0) 

7 - 2.0*3ETASQ*TRZSG  * 7 . Q*iETA *TRZ S3  - 9.C*TRZSQ) 

END  IF  (TYPE  OF  LOAD) 


5 G 0 C 
C 


RETURN 
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FUNCTION  OFT  NX  (TFiALX) 


C 

C 

C 

C 

C 

C 

C 

C 

r* 

U 

c 

c 


THIS  FUNCTION  CALCULATES  A VALUE  OF  THE  SEAM  ORIGINAL 
HINGE  LOCATION  EQUATION  FOR  THE  GIVEN  TRIAL  X (WHICH  IS  A 
OJESS  AT  THE  CORRECT  VALUE).  IF  THE  RESULT  RETURNED  IS 
ZERO,  THE  GUESS  IS  CORkECT  . 

THIS  FUNCTION  IS  CALLED  BY  THE  ROOT-FINDING  ROUTINE, 
BISECT.  IT  CALLS  NO  SUBFROC EDURES  OTHER  THAN  THE  LIBRARY 
FUNCTION,  EXF. 


INTEGER  BADXFG, 
COMMON  /FLAGS  / 

1 

COMMON  /INPUTS/ 


1 

2 

COMMON  /CONSTS/ 

1 

2 

3 

4 


BPFLAG,  DON  EFG , EOFLAG 

BADXFG,  BFFLAG,  DUNEFG , EOFlAG,  IEKRFG,  KOUNT, 

LUADFG,  MECFFGf  M1STFG,  NCONFG 

A,  AIDA,  ALPHA,  AR,  BETA,  D,  F,  H,  PC,  PE, 

Q , SIGMAC,  SIGMAR,  TAU,  T I NCR , f MAX,  TPRINT, 
XM,  WAVEFN 

ACU6E,  ARSQ,  ARSZP1,  ARZSP1,  ASQ,  B,  BARA1 , 
BARC1,  BARG1,  BARH1,  3ARU1,  BET  A Cb , BETAFR, 

BET  AS  G,  DELTAK,  EPSLNU,  EX  P BE  T , FOURTH,  HALF, 
GNEMZ , GNEPZ , SIXTH,  THETU1 , THIRD,  W,  XMU,  Z, 
ZCUbE,  ZFOUR,  ZSQ 


TRXSQ  = TRIA  LX*T  RI ALX 
ONE  M T X = 1.0  - TRIALX 

USE  EQUATION  DEFENDING  ON  TYPE  OF  LOAD  FOR  BEAM 
IF  LINEAR  LOAD 

IF  (LOAOFG. NE. 1 > OO  TO  10C0 

BFTNX  = AS  Q*T  RX  SQ*  (PE  - PC*ONEMTX  - AIDA*W)  - 6.G*F*XMU 


1~  BLAST  LOAD 

1000  IF  ( LOAUFG.NE. 2 ) GO  TO  5000 


BFTNX  = ONEMT  X*BETACB*  (ASUMRXSGMPE  - AIDA*W)  - 6.0*F*XMU) 

1 + 2.  G *ASQ*PCMEXF  < jE  i m*  ( -ONEMTX)  > 

2 *(3.0*ONEMTX*(BETASQ*TRXSQ  2.0 

3 - 2.0*6ETA*TRI  AlX) 

4 ♦ BETAMkXSQMBET  A*TRIALX  - 1.0) 

5 - 6. G*ONEMT X*EXP6ET 

a - BE T A*TRXSQ* ( 3ET A - 1.0)  ) 


END  OF  CASE  UN  TYPE  OF  LOAD 


50 OC  RETURN 
C 

END 


. 
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FUNCTION  FF  T NX  (TRIALX) 


w 

\ 


i 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


THIS  FUNCTION  CALCULATES  A VALUE  uF  THE  PLATE  ORIGINAL 
HINGE  LOCATION  EQUATION  FOR  THE  GIVEN  TRIAL  X (WHICH  IS  A 
GUESS  AT  THE  CORRECT  VALUE).  IF  THE  RESULT  RETURNED  IS 
ZERO,  THE  GUESS  IS  CORRECT. 


THIS  FUNCTION  IS  CALLED  BY  THE  ROOT -FINDING  PROCEDURE, 
BISECT.  IT  CALLS  THE  SU3PFCCEDURE,  CAL6AR,  TO  CALCULATE 
THE  ELEMENTS  OF  DELTA  DOT  DOT  AND  THETA  DOT  DOT  AT  TIME 
= 0,  FUR  THE  GIVEN  TRIAL  X. 


CUMMON  /BARS  / 

1 

2 

COMMON  /INPUTS/ 

1 

2 

COMMON  /CONSTS/ 

1 

2 

3 

'•* 


BARA,  BARB,  BARC,  BARD,  BARDEL,  BARDNM,  BARG, 
BARH,  BARJ,  BARK,  BARR,  BARS,  BART,  BARTHE,  BARU, 
EXPBX1,  EXP3ZX 

A,  AIDA,  ALPHA,  AR,  BETA,  D,  F,  H,  PC,  PE, 

Q,  SIGMAC,  SIGMAR,  TAU,  TI NCR , TMAX,  TPRINT, 

XM,  WAVEFN 

ACUBE,  ARSQ , ARSZP1,  ARZSP1,  ASQ,  3,  BARA1, 

BARCl,  BARG  1 , BARH1,  BAR J1 , BETAC3 , BETAFR, 

BE  T AS  C,  OELTAK,  EPSLNU,  EXP3ET , FOURTH,  HALF, 
ONEMZ , CNEPZ,  SIXTH,  THETU1,  THIRD,  W,  XMU,  Z, 
ZCU9E , ZFOUR , ZSQ 


CALL  CAL6AR( TRIALX) 

C 

PFTNX  = A * TRIALX  * (BARTHE  + B ARC*3  A RD  ♦ B ARG  *3  ARH ) 
1 - BARDEL  + DELTAK 

C 

RETURN 

C 

END 
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I 

i 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

0 

c 

c 

c 

c 


SUBKUJTINE  TSTEP  (LOOPFG) 

THIS  SUBROUTINt  CONTAINS  T HE  Li.uP  OF  OPERAi  IONS  PERF  ' EG 
0>J  THE  GIVEN  VARIABLES  FOR  THE  CASE  FO*  ONE  TIME  STEP. 

THIS  SUBROUTINE  IS  CALLEG  BY  Thi  T I ME -L  OOP- CD  NT  ROL 
PROCEDURE,  TCNTRL.  11  CALLS  T H:  FOLLOWING  SUbPROCEDOr-.ES : 
CHEKXHJ  TO  CHECK  THE  LOCATION  OF  THE  HINGE,  XH,  OUkING 
MECHANISM  TWO  CF  A CASE 

PRINTRI  TO  PRINT  OUT  A LINE  OF  RESULTS  FOR  THE  GIVEN 
TIME  STEP 

RUNGEKJ  TO  PERFORM  THE  RUNGE-KUTTA  COH  UTATIC  THE 

VARIABLES  AT  THE  CURRENT  TIME  STEF. 

THE  PARAMETER,  LOOPFG,  IS  THE  LOOP-CONT RDL-Fl AG , WHICH  IS  SET 
BY  THIS  ROUTINE  IF  THE  LOOF  SHOULD  BE  TERMINATEu  FUR  At,Y 
REASON. 

THE  CASE-IS -DONE-FLAG,  DONEFG,  CAN  BE  SET  IN  TWO  W AYS  t 

1.*  IF  EITHER  VELOCITY,  DELIA  DOT  OR  THETA  DOT,  GO  NEGATIVE 
INDICATING  THAT  THE  MAXIMUM  DEFLECTION  HAS  BEEN  REACHED 
DR  THE  PRESSURE  WAS  NCT  SUFFICENT  TO  SHOW  A RESPONSE 
AND  THE  CASE  CAN  BE  TERMINATED 
2 • t IF  THE  TIME  MAXIMUM  ( TMAX)  FOR  THE  CASE  TO  BE  RUN  HAS 
BEEN  REACHED,  IN  WHICH  CASE  THE  CASE  MUST  BE  ENDED. 


INTEGER  BAOXFG, 
COMMON  /FLAGS  / 

1 

COMMON  /INPUTS/ 

1 

COMMON  /RESULT/ 

1 

COMMON  /PRINTS/ 

1 

COMMON  /CONSTS/ 

1 

2 

3 

L 


5PFLAG , DONEFG,  EOFLAG 

3ADXFG,  BFFLAG,  UONEF  j , EOFlAG,  IERRFG,  KCUNT, 

LOALFG,  MECHFG,  M1STFG,  NCONFG 

A,  AIDA,  ALPHA,  AR,  BETA,  D,  F,  H,  PC,  PE, 

Q,  SIGMAC,  SIGMAR,  TAU,  T I NCR , TMAX,  T PRINT , 

XM,  WAVEFN 

DELTA,  T,  THETA,  THETAD,  V,  VEL,  WF,  WK,  WP , 

X , XH 

FLAG,  MAXLIN,  naida,  nf,  numlin,  numpag,  nwavef, 
STEPCT,  TIM.NOW,  T 1TLE  ( 15)  , TODAY,  TYPE 
ACUBE,  AF.SQ,  ArSZPI,  ARZSP1,  ASQ,  B,  BARAl, 
BARC1,  BARG  1 , 3AKH1,  3ARJ1  , BETACB,  BET  AFtt., 

BET  ASQ,  DELTAK,  EPSLNU,  EXP3ET,  FOURTH,  HALF, 
ONEMZ,  CNEFZ,  SIXTH,  THETU1,  THIRD,  W,  XMU , Z, 
ZCUBE , ZFCUH,  ZSQ 


DATA  STAR  /1H*/ 


COMPUTE  EQUATIONS  OF  MoTlUN  AND  WORK  FOR  THE  NEXT  TIME  STEP 

CALL  RUNGEK 
WK  =WF-WP-V 

T = T ♦ T1NC.R 

STcPCT  = STEPCT  + TINCI 


I ~ MECHANISM  1 

JF  ( MECHFG. NL. 1 ) GO  TO  1100 

16  6 . 


1 


c 

C THEN  COMPLETE  VEL  ANU  CELTA  CALCULATIONS 

C 

VEL  = THETAD  * A 
CELT  A = THETA  * A 


C 

C ELSE  CONTINUE 

c end  if  <mech  d 

c 

C IF  EITHER  VELOCITY  IS  ZERO 

C 

HOC  IF  TVfL.GE.0.0  .AND.  T HE  T Au  . G E . 0 . 0 ) GO  TO  1300 
C 

C THEN  WFUTE  MESSAGE  ANU  SET 

C CASE-I5-UCNE-FLAG  TO  TERMINATE  CASE 

C 

C IF  THIb  IS  THE  FIRST  TIME  STEP 

IF  (T.NE.TINCR)  GO  TO  1150 
C THEN  WRITE  NO  RESPONSE  MESSAGE 

WRITE  ( 6 , 9 9 0 3 ' 

GO  TO  1200 

C EL  WRITE  MAXIMUM  DEFLECTION  FOUND 

1 15  C WRITE  (6,9901)  DELTA,  T 

C 

C END  IF  (1ST  TIME  STEP) 

C 

1200  LOOPFG  = 1 
DONE  F G = 1 
WRITE  (6,9904) 

GO  TO  2000 
C 

C ELSE  CONTINUE  COMPUTATIONS 

C 

C IF  MECHANISM  2,  GET  NEW  HINGE  LOCATION  <XH) 

1300  IF  ( MECHFG .EQ. 2)  CALL  CHEKXH 
C 

C IF  HINGE  LOCATION  IS  OKAY 

IF  (6AOXFG.NE.0)  GO  TC  1500 


C 

C 

C 

C 


C 

C 

C 

C 

c 

c 


THEN  CHECK  FOR  FAILURE  WITH  THETA  U 

AND  PRINT  RESULTS  IF  END  OF  TIME  STEP  INTERVAL 

THETAU  = D * (SORT (THETU1*XH  * l.O)  - 1.0)  / XH 

IF  (THETA. GT. THETAU)  FLAG  = STAR 

IF  ( (STEPCT  T INCH)  .GT.TPRINT  ) CALL  PRINT  R 

IF  TIME  HAS  REACHED  THE  LIMIT  (TMAX ) 

IF  ( TMAX .GE. ( T + TINCR))  GO  TO  1400 

THEN  SET  LOOP -CONTROL-  AND  CASE-I S-DON E-Ft AGS 
AND  WkITE  TIME-tXCEEDEC  MESSAGE 
LOOPPS  = 1 
DONEPG  = 1 
WRITE  (6,9102) 


i 


I 
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WRITE  (b,  990*+) 

: else  continue 

END  IF  ( T =TM AX ) 

140C  GO  TO  1900 

ELSE  SET  LOOP  CONTROL  FLAG  TO  TERMINATE  THIS  TRY 
: BECAUSE  OF  A BAD  HINGE  LOCATION 

1500  LOOPFG  = 1 

END  IF  (GOOD  XH) 

1900  GO  TO  2000 

END  IF  (ZERO  VELOCITIES) 

2000  RETURN 

FORMAT  STATEMENTS 

9901  FORMAT  (1H0,  20X,  ’MAXIMUM  DEFLECTION  = ’,  G15.8,  * AT  TIME  = 

1 G15.3) 

9902  FORMAT  (1HO,  3 GX,  ’TIME  EXCEEDED’ ) 

9903  FORMAT  (1H0,  2GX,  ’INSUFFICIENT  PRESSURE  TO  GIVE  A RESPONSE’) 
9909  FORMAT  (1H0,  20X,  ’AN  ASTERISK  INDICATES  THAT  A REINFORCING  ’ 

1 ’ELEMENT  HAS  FRACTURED’) 


SUBROJT 1 NE  k U N G fc.  K 
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c 

c 

c 

n 

o 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 
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HIS  SUBROUTINE  SOLVES  FIRST  AND  SECOND  ORDER  SIMULTANEOUS 
DIFFERENTIAL  EQUATIONS  US  I N 0 THE  RUNGE-KJTTA  FOURTH-ORDER 
METHOD. 

THE  FIRST  ORDER  SOLUTIONS  OF  THE  DIFFERENTIAL  EQUATIONS  OF 
THE  HIGHEST  ORDER  (WHICH  CAN  BE  FIRST  OR  SECOND)  ARE  STORED 
IN  THE  ARRAY  CALLED  Y.  THE  FINAL  SOLUTIONS  OF  THE  DIFFER- 
ENTIAL EQUATIONS  (WHICH  ARE  ACTUALLY  THE  FIRST  ORDER  RESULTS 
STORED  IN  THE  ARRAY,  V)  Ar.E  STORED  IN  THE  A Rt  AY  U,  WITH 
MATCHING  SUBSCRIPT. 

THE  ARRAY,  ARG,  HOLDS  THE  VALUES  TO  BE  USED  IN  EACH  STEP 
OF  COMFUTING  T HE  FUNCTION.  THE  FIRST  ELEMENT  OF  ARG,  ARG(i), 
IS  ALWAYS  USED  FOR  THE  TIME  ARGUMENT,  T. 

THE  ARRAYS,  DX  AND  DXCX,  ARE  USED  FOR  THE  VALUES  OF  THE  FIRST 
AND  SECOND  ORDER  CHANGE  (DELTA)  DETERMINED  FOR  EACH  VARIABLE 
BY  THE  RJN  GE-KUTT  A METHOD.  THESE  ARE  ADDED  TO  EACH  ELEMENT 
OF  Y AND  U,  RESPECTIVELY,  TO  DETERMINE  THE  NEW  VALUES  OF 
THE  VARIABLES  AT  THE  END  OF  THE  TIME  S T E 3 . 

THIS  SUBROUTINE  CALLS  THE  SUBPROCEDURE,  COMP,  WHICH 
ACTUALLY  COMPUTES  THE  PROPER  FUNCTIONS  AT  EACH  POINT  OF  THE 
RUNGE-KUTTA  CALCULATION. 

THE  ARRAYS,  AC,  Al,  A2,  AND  A3,  HGlQ  THE  INTERMEDIATE 
RESULTS  OF  THE  RUNGE-KUTTA  CALCULATION,  SHOWING  THE  VALUE 
Or  EACH  FUNCTION  AT  EACH  POINT  OF  THE  CALCULATION.  THESE 
VALUES  ARE  COMBINED  TC  FORM  THE  VALUES  OF  THE  ARRAYS,  UX 
AND  DXDX , ACCORDING  TU  THE  RUNGE-KUTTA  FORMULA. 

FDR  FUKTHEk  INFORMATION  ON  THE  METHOD  USED  HERE,  THE  READER 
IS  DIRECTED  TO  THE  TEXT,  INTRODUCTION  TO  NUMERICAL  ANALYSIS 
BY  HILDEBRAND  (NEW  YuRKS  MCGRAW-HILL,  1956),  FP.  233-239. 

THIS  SUBROUTINE  IS  CALLED  BY  ThE  PROCEDURE,  TSTEP.  IT 
CALLS  THE  SUGFRCCEDURE,  CCMP. 


INTEGER  BAD  XFG , 
CUMMON  /FLAoS  / 

1 

COMMON  /INPUTS/ 

1 

7 

COMMON  /RESULT/ 

1 


BFFLAG,  DON  EFG , cuFLAG 

BADXFG,  BFFLAG,  DGNEFu,  EOFlAG,  IERRFG,  KOUNT, 

LOAUFG,  MEC  Hr G,  MlSTFG,  NCONFG 

A,  AIDA,  ALPHA,  AR,  BETA,  D,  F,  H,  PC,  PE, 

0,  SIGMAC,  SIGMAR,  TAU,  T1NCR,  TMAX,  TPkINT, 

X M , WAVEFN 

DELTA,  T,  THETA,  THETAO,  V,  VEL,  WF,  WK,  WF, 

X , XH 


DIMENSION  Y ( 5 ) , U ( 2 ) , 

1 D X ( 5 ) , 0X0 X ( 2 ) , 

2 AC(5),  Al(5>,  A2 ( 5 ) , A3(6), 

j A RG  ( 5 ) 
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INITIALIZE  VALUES  OF  FIRST  AND  SECOND  ORDER  EQUATION 

Y ( 1 ) = VEL 
U(l)  = DELTA 

Y ( 2 ) = THETAD 
U ( 2)  = THETA 
Y ( 3)  = WF 

Y ( 4 ) = WP 
Y (E)  = V 

SET  UP  ARGUMENTS  FOR  FIRST  STEP  OF  RUNGE-KUTTA 

ARG(l)  = T 
ARG ( 2)  = DELTA 

AKG ( 3 ) = VEL 
ARG(4>  = THETA 
ARG(i>)  = THETAD 
CALL  COMP  (ARG,  A 0 > 

SET  UP  ARGUMENTS  FOR  THE  SECOND  STEP 

HALFTI  = 0 • 5 *T I NCR 
ARG  ( 1)  = T + HALFTi 

ARG ( 2 > = DELTA  + HALFTI*VEL 
ARG ( 3>  = VEL  + 0. 5*AG (1) 

Ai\GU)  = THETA  *■  HAL  FT 1*T  HET  AH 
A F G ( 5 1 = THETAD  ♦ Q.5*AQ(2) 

CALL  COMP  (ARG,  Al) 

SET  UP  ARGUMENTS  FOR  THE  THIRD  STEP 

ARG ( 2)  = DELTA  * H ALFT I* V IL  + 0 . 5* H AL FTI * AD ( 1 ) 

ARG  ( 3)  = VEL  «•  0.5*A1(1) 

A RG ( 4)  = ThETA  + H AL f T I* T HET AD  + G . 5*HALFTI* AO ( 2) 

ARG  ( 5)  = THETAO  + i).5*Al(2) 

CALL  C CM  P (ARG,  A2) 

SET  UP  ARGUMENTS  FOR  THE  FOURTH  AND  LAST  STEP 
ARG ( 1)  = T * T I NCR 

ARG  ( 2 ) = DELTA  4-  T 1 NCR*  V EL  4 H A L FI  I * A 1 ( 1 > 

AKG(  3)  = VEL  ♦ A 2 ( 1 ) 

ARG  ( 4)  = THETA  ♦ T INCR*T  HET  AD  + HA L F T I * A 1 ( 2 ) 

ARG ( 5 ) = THETAD  f A2(2) 

CALL  COMP  (ARG,  A3) 

PUT  PIECES  TOGETHER 

DO  10J0  1=1,2 

UXDX(I>  = TINCR*Y(I)  + TINCRMAG  (I)  *A  1 ( I > ♦ A 2 ( I ) )/6.0 
1GOC  U(I)  = U (I)  + DXUX  (I) 

DO  11JU  I = l,a 


J 


ux  C X > = < AO ( I ) ♦ 2.0*A1(I)  + 2 • 0 * A 2 ( I ) ♦ A3(I))/6.0 

y (i)  = y ( I ) ♦ oxu > 

COMPUTATION  COMPLETE 

\l  EL  = Y (1) 

DELTA  = U<1> 

THETAD  = Y ( 2 ) 

THETA  = U ( 2) 

WF  = Y (3) 

HP  = Y ( 4) 

\l  = Y ( 5 ) 

RETURN 

END 
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C 

C 

C 
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C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 
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SUbROJT INE  COMP  (ARG,  AA) 


TH 

IS 

SU3ROUTI 

SE 

TS 

OF 

FUNCT 

EN 

T I 

al 

EQUATI 

S£ 

T 

OF 

EQUATI 

TY 

PE 

OF 

LOAD 

0- 

T 

HE 

CASE. 

NE  ACTUALLY  COMPUTES  ONE  OF  EIGHT  POSSIBLE 
IONS  WHICH  DESCRIBE  THE  SIMULTANEOUS  DIFFER- 
C NS  BEING  SOLVED  FOR  EACH  TIME  STEP.  THE 
ONS  USED  DEFENDS  ON  THE  TYPE  OF  CASE  AND 
AS  WELL  AS  THE  CURRENT  STATE  (OR  MECHANISM) 


THE  INPUT  AND  OUTPUT  PARAMETER  ARRAYS  ARE  DEFINED  AS 
FOLLOWS* 

ARG*  INPUT  TO  THE  FRCCECURE,  CONTAINS  THE  ARGUMENTS 
BE  USED  IN  COMPUTING  THE  FUNCTIONS 
A A I OUTPUT  RESULTS  TO  BE  RETURNED  TO  THE  CALLING 

PROCEDURE,  RUNGEK. 


TO 


THIS  SUBROUTINE  CALLS  THE  SUBPROCEDURE,  CAlBAR,  TO  COMPUTE 
PORTIONS  OF  THE  PLATE  EQUATIONS  OF  MOTION.  IT  ALSO  USES  THE 
SYSTEM  LIBRARY  ROUTINE,  EXP.  THIS  ROUTINE  IS  CALLED  BY  THE 
RJNGE-KUTTA  SUBPROCEOUF  E , RUNGEK. 


1 

1 

2 

1 

2 

1 

2 

1 

1 

2 

3 

k 


INTEGER 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 


BAOXFG  , 

/flags  / 
/BARS  / 

/COMPS  / 


/INFUTS/ 


/RESULT/ 

/CONSTS/ 


BPFLAG,  DONEFG,  EOFLAG 

BADXFG , BFFLAG,  DONEFG,  EOFLAG,  IERRFG,  KOUNT , 
LOAOFG*  MECHFG,  MlSTFG,  NCONFG 
BARA,  BARB,  BARC,  BARD,  BARDEL,  BARDNM , BARG, 
BARH,  3ARJ,  BARK,  BARR,  BARS,  BART,  BARTHE,  BARU, 
EXFQX1,  EXFBZX 

ATHEDB,  ATHEDL , ATHED1,  ATHED2,  AVDOT,  AWFDB, 
AWFUL,  BTHEB2,  BTHED3,  BVDO  T , BWFDOT , BWPDOT, 
PWFDB , PWFOL,  PWPDOT,  PVDOT 

A,  AIDA,  ALPHA,  AR,  BETA,  D,  F,  H,  PC,  PE, 

Q,  SIGMAC,  SIGMAR,  TAU,  TINDR,  TMAX,  TPRINT , 

XM,  WAVEFN 

DELTA,  T,  THETA,  THETAD,  V,  VEL,  WF,  WK,  WP, 

X,  XH 

ACU3E  , ARSQ , ARSZP1,  ARZSP1,  ASQ,  B,  BAkAl, 

BARC1,  BARG  1 , BARH1,  BARJ1 , 3ETACB,  BETAFR, 

B ET ASQ,  OELTAK,  EFSLNU,  EXP3ET  , FOURTH,  HALF, 
ONEMZ  , CNEPZ,  SIXTH,  THETU1,  THIRO,  W,  XMU,  Z, 

Z CUBE , ZFOUR , ZSQ 


DIMENSION  A RG ( 5 ) , 


AA  (5) 


DETERM  IN 
USING 


E FTNT 

CASE  ON  RATIO 


OF 


T * TAU 


WITH  THE  WAVE-FUNCTION 


FTNT 
r RAT 
IF  ( 

1 

IF  ( 

1 


= 0.0 


ID  = ARG ( 1 ) / 

TAU 

TRATIO. LE. 1.0 

.AND  . 

WAVEFN. EG. 

r TN  T = (1.0 

- TRATI 

C)  * EX P (- 

TRATIO. LE. 1.0 

.AND. 

WAVEFN. EQ. 

FTNT  = 1.0 

1.0  ) 
ALPHA 
2.0) 


* TRATIO) 


oo  ooo  ooooo  ooooooo  ooooooooo  oooo 


IF  FIRST  TIME  OR  SECOND  TIME  THROUGH  THIS  ROUTINE 


IF  (H1STFG.EQ.0)  GO  TO  900 

THEN,  TO  AVOID  ZERO  DIVISION  AND  TO  INITIALIZE  THE  FIRST 
TIME  STEP  FOR  THE  RUNGE-KUTTA  METHOD,  SET  THE  SECOND  AND 
FOURTH  ARGUMENTS  SUCH  THAT  CUKRX  WILL  REFLECT  THE  VALUE 
OF  THE  ORIGINAL  HINGE  LOCATION.  THIS  WILL  BE  DONE  THE 
FIRST  AND  SECOND  TIME  THE  PROCEDURE  CGMP  IS  CALLED  » 

THAT  IS,  FOR  THE  CALCULATION  OF  THE  FIRST  TWO  PARAMETERS 
OF  1 HE  RUNGE-KUTTA  FOR  THE  FIRST  TIME  STEP  ONLY. 

ARG( 2 ) = XH 
ARG( 4)  = 1.0 
M1STFG  = M 1ST  FG  - 1 

ELSE  CONTINUE 

END  IF  (1ST  OR  SECOND  TIME  ThRU) 


CASE  ON  COMBINATION  OF  LOAD-,  MECHANISM-,  AND  BE AM- PL AT E -FL AGS 

9 0 C LMBPSN  = (LOADFG  - 1)*4  * (BPFLAG  - 1>*2  + (MECHFG  - 1)*1  «■  1 
GO  TO  (1000,  2000,  3000,  4000, 

1 5000,  6000,  7000,  8000),  LMBPSM 

LINEAR  LOAD  CASES 

THIS  IS  A BEAM  CASE  IN  MECHANISM  1,  LINEAR  LOAD 

1000  CURRX  = 1.0 
A A ( 1 ) = 0.0 

AA(2>  = FT  NT  * ATHEDL  * ATHED2 
A A ( 3 ) = FT  NT  * ARG ( 5 ) * AWFDL 
AA  (4)  = BWPDOT  * ARG ( 5 ) 

A A ( 5 ) = AVDOT  * ARG ( 5 ) 

GO  TO  9000 

THIS  IS  A BEAM  CASE  IN  MECHANISM  2,  LINEAR  LOAD 

2000  CURRX  = ARG( 2 ) / (A  * ARG ( 4 ) ) 

A A ( 1 ) = FT  NT*  (PE  + FC* HALF* ( CURRX  * 1.0)  >/XM  - DELTAK 
AA  (2 ) = FT  NT  * (3.0*PE  * 2. 0 *CURRX*PC) 

1 / (2.0  * A * CURRX  * XM) 

2 - BT HED2/CURRX  - 3THED3/ ( CURRX**3) 

A A ( 3 ) = F T NT  * ARG ( 5 ) * BWFDCT  * CURRX 

1 * (PE  * (1.0  - CURRX* HALF) 

2 ♦ PC  * HALF  * (1.0  - THIRD* CURRX* * 2 ) ) 

A A ( 4 ) = BWPDOT  * ARG  ( 5 ) 

A A ( 5 ) = cVDOT  * ARG( 5)  * CURRX  * (1.0  - CURRX*HALF> 

GO  ro  9000 

THIS  IS  A PLATE  CASE  IN  MECHANISM  1,  LINEAR  LOAD 
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3 0 0 C CURRX  = 1.0 

A A ( 1 ) = 0.0 

AA<2)  = FT NT  * ATHEOl  ♦ A THE C2 
A A ( 3 ) = F T NT  * A k G ( 5 ) * AKFOL 
AA (4 ) = F N POO  T * AkG ( 5 ) 

AA(5)  = A\/  DOT  * ARG  (5 } 

go  ro  gooo 

C 

C THIS  IS  A PLATE  CASE  IN  MECHANISM  2,  LINEAR  LOAD 

C 

4000  CURRX  = ARG ( 2 ) / (A  * ARG(4)> 

CALL  CALbA  R( CURRX ) 

AA (1 ) = FT  NT  * EARDEL  - DELTAK 

A A ( 2 ) = FT  NT  * BARTHE  BARC  * BARO  ♦ BARG  * BARH 

A A ( 3 ) = FT  NT  * ARG ( 5 ) * PWFDL 

1 * (PE*CURRX* (1 .0  - C'JRRX*HALF*ONEPZ  ♦ Z*THIRD*CURRX**2) 

2 + PC*CURRX* (THIRD  - ON EPZ*S IXT H * CU RRX * CURRX 

3 + (SIXTH  - Z*QN£MZ * FOU RT H ) *CURRX**3)  ) 

AA  (4  ) = PWFDOT  * ARG  ( 5 ) 

A A (5 ) = PO  DOT  * ARG  < 5)  * CURRX 

* (1.0  - ONEPZ*HALF*CURRX  ♦ Z*TH I RD * CU RRX*CURRX ) 

GO  TO  9000 

BLAST  LOAD  CASES 

THIS  IS  A BEAM  CASE  IN  MECHANISM  1,  BLAST  LOAD 

CURRX  =1.0 

A A ( 1 ) = 0.0 

A A ( 2 ) = F T NT  * ATHEDB  - ATHED2 
A A ( 3 ) = FT  NT  * ARG  ( 5 ) * AWFDB 
A A ( 4 ) = BWPDOT  * ARG ( p ) 

A A (5 ) = AtfDOT  * ARG( 5) 

GO  TO  9000 

THIS  IS  A BEAM  CASE  IN  MECHANISM  2,  BLAST  LOAD 

CURRX  = ARG ( 2 ) / (A  * ARG( 4 ) ) 

BETAX  = BETA  * CURRX 
ONEMCX  = 1.0  - CURRX 
EXP3CX  = EXP  (BETA  * (-ONEMCX)  ) 

A A ( 1 ) = ( F TNT /KM)  * (PE 

♦ (PC/ ( BET ASQ*GNEMCX) ) 

♦(BETA  - 1.0  - (BETAX  - 1.0)*EXPBCX)  ) 

- DELTAK 

AA ( 2 ) = ( 3 • 0 ♦ FTNT  / (XM* A*CUKRX**3 ) ) 

* (PE*HALF*CURRX*CUkRX 

♦ (PC/BETACB)  ♦ (EXPBCX*  (BETAXMBETAX  - 2.0)  ♦ 2.0) 

- 2. 0 ♦EXPBET)  ) 

- BTHED2/CURRX  - BTHEB3/ ( CURRX **3 ) 

A A ( 3 ) = FTNT  ♦ AP.G(5)  ♦ 3WFDGT 

♦ (PE  ♦ CURRX  *(1.0  - CURRX*HALF) 

♦ ( PC/OET  AC6)  * (EXFBCX*  (2. 0 - BETAX) 

♦ BET  AX* (BET  A - 1.0)  - 2.0*EXPBET)  ) 


oooooo  ooo  ooo 


I 

( 


7 0 0 C 


80  00 


1 

2 

3 

4 

5 

6 

7 

8 


A A ( 4 ) 

BWPDOT 

* ARG (5) 

AA  (5  ) 

= 

DVDOT  * 

ARG < 5 ) 

* CURRX  * (1. 

0 

go  ro 

9000 

THIS 

IS  A PL 

ATE  CASE 

IN  MECHANISM 

1 

CURRX 

s: 

1.  0 

A A ( 1 ) 

is 

0.  0 

AA  (2 ) 

=: 

FT  NT  * 

ATHEO  1 

♦ ATHED2 

AA  (3  ) 

s: 

FT  NT  * 

ARG  (5 ) * 

AWFDB 

A A ( 4 ) 

s 

FH  FOOT 

* AkG ( 5 ) 

AA  ( 5 ) 

s 

Al/OOT  * 

ARG (5) 

GO  TO 

9000 

THIS 

IS  A FLATE  CASE 

IN  MECHANISM 

2 

CURRX 

= 

A R G ( 2 ) 

/ (A  * ARG ( 4 ) ) 

CAL. 

CALBAR(CJRRX) 

AA  (1) 

FTNT  * 

eAROEL 

- DELTAK 

AA  ( 2 ) 

= 

FTNT  * 

BARTHE 

♦ BARC  * BARO 

AA  (3  ) 

FTNT  * 

ARG  (5  ) * 

PWFDB 

CURRX*HALF) 


BLAST  L 0 AO 


BARG 


BARH 


A A ( 4 ) 
A A ( 5 ) 


* (PE*BETAFR*ZSG¥  (CU«.RX  - Z*CURRX*CJRRX*HALF 

+ 0NEPZ*THIRD*CURRX**3) 

♦ PCM2.0*BETASQ'»-ZSQ,f-CURRX  - 4. 0 *3ETA*ZSQ*CURRX 

♦ EXP3ZX* (BETASQ*ZSQ*CURRX* (CURRX  - 1.0) 

- 3.0*BETA*Z*CURRX  «•  2. 0*BETA*Z  +■  3.0) 

♦ EXPBX1*ZSQMBETASQ*CURRX*C'JRRX  - 5 . 0 * BET  A*CURRX 

- BET ASQ*CURRX  + 2.0*BETA  ♦ 9.0) 

- EXP3ET*(2.0*BETA*Z*ONEPZ  ♦ 9.0*ZSQ  +■  3.0)  ) ) 

= PM  FOOT  * ARG ( 5 ) 

= Pi/  DOT  * ARG(  5)  ♦ CURRX 

* (1.0  - ONEPZ*HALF*CURRX  ♦ Z*  THIRD*  CURRX*  CURRX  ) 


END  CASE  ON  LOADFG,  BPFLAG  g,  MECHFG 


MULTIPLY  EACH  TERM  BY  THE  TIME  STEP  INCREMENT 

90 0 G DO  9100  1=1,5 

9100  A A ( I ) = TIMOR  * AA(I) 

C 

RETURN 

END 
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SUBROUTINE  CAtBAK  (XX) 


C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 


TH 

I 

S SUB 

ROUTINE 

CAL  CUL 

ATES 

PIE 

OF 

MOT  10 

N, 

THE 

BARRED 

ELEM 

ENTS 

THE 

TA  DO 

T 

DOT. 

TH 

I 

S PRO 

CE 

DURE 

IS 

CALL 

ED  B 

Y TH 

PTZERO 

• 

• 

TO  G 

£T 

ELEM 

ENTS 

OF 

;l  of  the 

MAKING  UP 


PLATE  EQUATIONS 
DELTA  DOT  DOT  AND 


MECHANISM  1 


£ FOLLOWING  SUBPROCEDURES: 
THE  EQUATIONS  OF  MOTION  TO 
CONSTANTS  FOR  THE  PROCEDURE, 


COMP: 


PFTNX: 


OF 

THE 

EQUATIONS 

OF 

MOTION 

TEP 

FOR 

THE  RUNGE- 

KUTTA 

OF 

THE 

E QUAT IONS 

CF 

MOTION 

NAL 

H I KG 

E LOCATION 

FUNCTION 

DEFINE 
COMP 

TO  COMPUTE  ELEMENTS 
AT  THE  GIVEN  TIME  S 
ROUTINE,  RUNGEK 
TO  COMPUTE  ELEMENTS 
MAKING  UP  THE  ORIGI 
AT  TIME  ZERO  IN  ORDER  TO  FIND  A ROOT  FOR  SAID 
FUNCTION, 

IT  CALLS  NO  SUBPROCEDURE  OTHER  THAN  THE  LIBRARY  SYSTEM 
EXPONENTIAL  FUNCTION,  EXP, 

ALL  INPUT  COMES  THROUGH  THE  COMMON  BLOCK,  CONSTS,  WITH 
THE  EXCEPTION  OF  THE  SINGLE  PARAMETER,  XX,  WHICH  IS  DEFINED 
AS  FOLLOWS: 

i:  THE  CURRENT  VALUE  CF  THE  HINGE  LOCATION  (FROM  COMP) 

2:  THE  CURRENT  GUESS  AT  THE  ORIGINAL  HINGE  LOCATION 
(FROM  PFTNX) 

3t  THE  FINAL  VALUE  OF  THE  HINGE  LOCATION  (FROM  PTZERO). 
Al-L  OUTPUT  IS  PASSED  THROUGH  THE  COMMON  BLOCK,  OARS. 


INTEGER  BADXFG,  BPFLAG,  OONEFG,  EOFLAG 
COMMON  /FLAGS  / BADXFG,  BPFLAG,  DONEFG, 

LOADFG,  MECHFG,  M1STFG, 
BARA,  BARB,  BARC,  BARB, 
BARH,  BARJ,  BARK,  BARR, 


COMMON  /BARS 


EOFLAG,  IERRFG,  KOUNT, 
NCCNFG 

BARBEL , BARDNM,  BARG, 
BARS,  BART,  BARTHE,  BARU, 


COMMON  /CONSTS/ 


COMMON  /INFUTS/ 


EXPBX1, 

EXP3 

ACUBE, 

ARSQ , 

BARC1, 

BARG1 

BET  AS  Q, 

DELT 

ONEMZ  , 

ONEPZ 

Z CUBE  , 

ZF  OUR 

A,  AIDA 

, ALP 

AKSZP1,  ARZSP1,  ASQ,  B,  BARA1, 
3ARH1,  BARJi,  BETACB,  BETAFR, 

, K,  EFSLNU,  EXP  BE  T , FOURTH,  HALF, 
SIXTH,  THETU1,  THIRD,  W,  XMU,  Z, 
ZSQ 


Q, 

XM  , 


SIGMAC, 

WAVEFN 


, AR,  BETA,  D,  F, 
SIGMAR,  TAU,  T I NCR , 


H,  PC,  PE, 
f MAX , TPRINT 


XXSQ  = 

XX 

*xx 

XXCUBE 

XX*  * 3 

ONEMXX 

= 

1.0  - XX 

onemz< 

s 

1.0  - z*x 

EXPBX1 

5 

EXP (BETA* 

EXPBZX 

- 

EXP  (BETA* 

3 ARD.NN 

= 

ARZSP1 * ( r 

= ARZSP1M  THIRD  - XX*FOUKTH)  ♦ ONE MZ*XX* FOURTH 


17b  . 


ooooo  oooooo  oooooo  oooooo 


C 


3ARC  = 

BARC1 

/ 

XX 

BARU  = 

( ARZSPl 

* (HALF 

BARG  = 

6ARG1 

/ 

XX  CUBE 

BAkH  = 

BARH1 

/ 

BARDNM 

XX*THIRD)  +•  ONEMZ*XX*THIRD ) / BARDNM 


CASE  ON  TYPE  OF  LOAD 
LINEAR  LOAD 

IF  (LOADFG.NE. i)  GO  TO  2000 

BARR  = 1.0  / (XM  * A * XX  * BARDNM) 

BARS  = PE*  (ARZSF1* (HALF  - XX*THIRD)  ♦ ONEMZ*  XX*T  NIRO) 

1 «•  PC*(ZCUBE*XX*AR*(THIRO  - XX*F0URTH) 

2 «■  XXSQ*  FOURTH* (ll.G*SIXTH  - ZSQ*  HALF  - XX)  ) 

BARTHE  = BARR  * BARS 

IF  HINGE  NOT  IN  CENTER  OF  PLATE 

IF  ( XX. EQ. 1.0)  GO  TO  2000 

THEN  COMPUTE  DELTA  DOT  DOT  ELEMENTS 

BART  = 1.0  / (XM  * ONEMZ X ) 

8ARU  = PE*  O NE  MZ X ♦ PC* ( THIRD* ( 1 . 0 * XX) 

1 - XXSQ*HALF* (ZSQ  * THIRD)  ) 

8ARDEL  = BART  * BARU 

ELSE  CONTINUE 
END  IF  (HINGE  IN  CENTER) 

IF  BLAST  LOAD 

2000  IF  (LOADFG.NE. 2)  GO  TO  5000 

BARA  = BAR A1  / BARDNM 
BAR 3 = PE*BETAFR*ZSQ*XXSQ* 

( AR*ZSQ*  ( HALF  - XX*THIRD)  ♦ ONEM  Z*  THI RD*  XX  1.0) 

♦ PC*(-EXPBET* (2.0*BETA*ZSQ*(AR  + 1,3)  ♦ 9.0*ZSQ 

♦ 6.0* AR* L - 3.0) 

* EXPBX1* (BETACB*Z3Q*XXSQ*0NEMXX 

♦ BET ASQ*ZSQ* (4.0*XXSQ  - 2.0*XX> 

- 9.0*3ETA*ZSQ*XX  + 2.0*d£TA*ZSU  «■  9.0*ZSC) 

♦ EXPBZX* (AR* (BETACB*ZFCUR*XXSQ*ONEMXX 
* BETASQ*ZCUBE* (3. C*XXSQ  - 2.0*XX> 

- 6.0*BET  A*ZSQ*XX  * 2.0*BETA*ZSQ  + 6.0*Z) 

- BETASG*ZSU*XXSQ 
♦ 3.0  *BET  A*  Z*XX  - 3.0)  ) 

BARTHE  = BARA  * BARB 

IF  HINGE  NOT  IN  THE  CENTER  OF  THE  PLATE 
IF  (XX. EQ. 1.0)  GO  TO  5000 

THEN  COMPUTE  DELTA  OCT  OUT  ELEMENTS 


1 

2 

3 

4 

5 
5 

7 

8 
9 
A 
3 


177 


c 

c 

c 

c 

c 

c 

c 


BARJ  = 3ARJ1  / (ONF.MXX  * 
BARK  = SET  ACB*PE*(Z*XXSQ 
♦ PC *(2.0* BETA  - 


BARDEL  = 


ONEMZ  x ) 

- ONEPZ*XX  ♦ 1,3) 

4.0 

♦ EXP8ZX*BETA*0NEMXXM 1.0  - BETA*Z*XX) 

♦ EXPBX1M4.0  + SETA  - 3.0*BETA*XX 

- BETASQ*XX*ONEMXX)  ) 

BARJ  * BARK 


ELSE  CONTINUE 
END  IF  (HINGE  IN  CENTER) 

END  CASE  ON  TYPE  OF  LOAD 


5000 


RETURN 

END 
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SUURUJT INE  CHEKXH 
C 

C HIS  SUBROUTINE  CHECKS  THE  VALUE  OF  THE  HINGE  LOCATION,  XH, 

C TO  SEE  THAT  IT  IS  WITHIN  RANGE,  AND  TO  SEE  IF  IT  HAS  MOVE U 

C FROM  A MECHANISM  2 To  A ML  IAN1SM  1 POSITION. 

C 

C IF  THE  HINGE  LOCATION  IS  NEGATIVE  (OFF  THE  END  OF  THE  BEAM  OR 

c plate),  the  data  must  be  bad.  if  the  hinge  location  has  moved 

C WITHIN  TWO  PERCENT  OF  THE  CENTER  OF  THE  3EAM  OR  PLATE,  IT  IS 

C CONSIDERED  TO  BE  AT  THE  CENTER,  WHICH  IS  THE  MECHANISM  1 

C POSITIUN,  AND  THE  MECHANISM  FLAG  IS  CHANGED  TO  REFLECT  THIS, 

C IF  THE  HINGE  LOCATION  HAS  MOVED  FAST  THE  CENTER  OF  THE  BEAM 

C OR  PLATE  BY  MGRE  THAN  TWO  PERCENT,  THE  COMPUTATION  IS  CON- 

C SIDERED  TOO  INACCURATE,  AND  SO  THE  TIME  INCREMENT,  T I NCR , 

C IS  HALVED,  A MESSAGE  IS  PRINTED,  AND  FLAGS  ARE  SET  TO  SHOW 

C THAT  THIS  HAS  HAPPENED  AND  TO  RUN  THROUGH  THE  COMPUTATION 

C AGAIN  WITH  THE  SMALLER  TIME  STEP.  THE  MAXIMUM  NUMBER  OF 

C TIMES  WHICH  THIS  CAN  BE  ATTEMPTED  IS  TWO  (AND  IS  STORED  IN 

C THE  CONSTANT,  MAXTRI,  IN  THE  SUBROUTINE,  TCNTRL). 

C 

C THE  BAD-HINGE-LOCATICN-FLAG,  BADXFG,  IS  SET  IF  THE  HINGE 

C LOCATION  FOUND  IS  NEGATIVE  OR  PAST  THE  CENTER  POINT  BY 

C MORE  THAN  TWO  PERCENT. 

C 

C THE  CASE -IS -DONE- FLAG,  OQNEFG,  IS  SET  IF  THE  HINGE  LOCATION 

C IS  NEGATIVE. 

C 

C THE  MECHANISM-FLAG,  MECHFG , IS  CHANGED  FROM  2 TO  1 IF  THE 

C HINGE  LOCATION  HAS  MOVED  TO  THE  CENTER  OF  THE  BEAM  OF  PLATE. 

C 

C THIS  SUBROUTINE  IS  CALLED  EY  THE  FOLLOWING  PROCEDURES! 

C TCNTRL!  TO  CHECK  THE  ORIGINAL  HINGE  LOCATION  AT  T=0 

C TSTEP t TO  CONFUTE  AND  CHECK  THE  HINGE  LOCATION  AT 

C EACH  SUCCESSIVE  TIME  STEP. 

C THIS  RUU  TINE  CALLS  NO  SUBPROCEDURES. 

C 

C 

BPFLAG,  DONEFG,  LGFlAG 

BADXFG,  BPFLAG,  DONEFG,  EOFLAG,  IERRFG,  KOUN I , 
LOADFG,  MECHFG,  MiSTFG,  NCONrG 
A,  AIDA,  ALPHA,  Ak,  BETA,  D,  F,  H,  PC,  PE. 

Q,  SIGMAC,  SIGMAR,  TAU,  TINCR,  TMAX,  TPRInT  , 

XM,  WAVEFN 

DELTA,  T,  THETA,  THETAO,  V,  VEL,  WF,  WK,  * , 

X,  XH 


INTEGER  BADXFG, 
COMMON  /FLAGS  / 

1 

COMMON  /INPUTS/ 

1 

2 

COMMON  /RESULT/ 

1 


IF  THIS  IS  NOT  THE  FIRST  TIME  STEP  (T/=0) 


k 


IF  (T.EQ.0.0  .OR.  THETA. EQ.O. 0)  GO  TO  10U 

THEN  WE  NEED  TC  COMFUTE  HINGE  LOCATION  FROM 
XH  = DELTA  / THETA 
X = XH  / A 
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MICROCOPY  RESOLUTION  TEST  CHA.RT 

NAIlONAl  BURLAU  Of  STANDARDS  1963-^ 


C IF  THE  HINGE  LOCATION  HAS  NOT  REACHED  THE  FINAL  HINGE  LOC 

100  IF  (X.GT.0,98)  GO  TO  3000 
C 

C THEN  CASE  IS  STILL  WITHIN  MECHANISM  2 

C IF  HINGE  LOCATION  IS  NEGATIVE 

C 

IF  (XH.GT. 0. 0)  GO  TO  2000 
C 

C THEN  DATA  MUST  BE  BAD  * WRITE  ERROR  MESSAGE  AND 

C SET  BAO-HINGE-  AND  C ASE-IS-DONE-  FLAGS  TO  TERMINATE  CASE 

C 

WRITE  (6,9901) 

BAOXFG  = 1 
DONEFG  = 1 
C 

C ELSE  CONTINUE — HINGE  LOCATION  IS  WITHIN  MECH  2 

C 

C END  IF  (NEG  XH) 

C 

2000  GO  TO  5000 
C 

C ELSE  HINGE  LOCATION  HAS  MOVED  TO  MECHANISM  1 LOCATION 

C NEED  TO  CHECK  RANGE 

C 

C - IF  HINGE  LOCATION  IS  WITHIN  2 7.  OF  THE  FINAL  HINGE  LOCATION 

3000  IF  (X.GT.1.02)  GO  TO  4000 
C 

C THEN  SET  HINGE  TO  FINAL  HINGE  LOC  AND  SET  FLAGS  TO  SHOW 

C MECHANISM  1 HAS  BEEN  SUCCESSFULLY  REACHED 

C 

XH  = A 
X = 1.0 
MECHFG  = 1 
GO  TO  4999 
C 

C ELSE  HINGE  LOCATION  HAS  OVERSHOT  THE  FINAL  HINGE  LOC 

C NEED  TO  DECREASE  TIME  STEP  ANO  TRY  AGAIN  FROM  T=0 

C 

4000  TINCR  = TINCR/2.0 
fcADXFG  = 1 

WRITE  (6,9900)  A,  XH 
C 

C END  IF  (27.  OF  FINAL  HINGE  LOC) 

C 

4999  GO  TO  5000 
C 

C END  IF  (XH  REACHED  FINAL  HINGE  LOC) 

C 

5000  RETURN 
C 

C FORMAT  STATEMENTS 

C 
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9900  FORM  AT  (1H0,  *HINGE  LOCATION  HAS  OVERSHOT  FINAL  HINGE  LGC*/ 

1 * FINAL  HINGE  LOC  = *,  G15.8,  * HINGE  IS  AT  *,  G15. 8/ 

2 * TIME  INCREMENT  HAS  BEEN  HALVEO— CASE  WILL  BE  RERUN*) 

9901  FORMAT  C1H0,  * HI NGE  LOCATION  IS  NEGATIVE- -CASE  IS  TERMINATED*/ 

1 * CHECK  INPUT  DATA  VALUES*) 

END 

* 

3 


ooo  o ooooooo  oooooo  oooooooo 


SUBROUTINE  PRINTR 


THIS  SUBROUTINE  PRINTS  THE  LINE  OF  RESULTS  FOR  THE  TIME,  T, 
IT  IS  CALLED  BY  THE  SUBPROCE CUKES , T2EK0  ANO  TSTEP.  IT 
CALLS  SUBPROCECURE  PAGE  TO  HEAD  A NEW  PAGE,  IF  THE 
CURRENT  LINE  COUNT  FOR  THE  PRINTEO  PAGE  EXCEEDS  THE 
MAXIMUM  NUMBER  OF  LINES  PER  PAGE. 


COMMON  /PRINTS/ 

1 

COMMON  /RESULT/ 

1 


FLAG,  MAXLIN,  NAIDA,  NF,  NUNLIN,  NUMPAG , NWAVEF, 
STEPCT,  TIMNOW,  TITLE115),  TODAY,  TYPE 
DELTA,  T,  THETA,  THETAD,  V,  VEL,  WF,  WK,  WP, 

X,  XH 


IF  LINE  COUNT  EXCEEDS  MAX/FAGE 

THEN  HE  AO  A NEW  PAGE 

IF  < NJMLIN.3E. MAXLIN)  CALL  PAGE 

ELSE  CONTINUE 
END  IF  (LINECOUNT) 

P<INT  LINE  OF  RESULTS,  AOD  1 TO  LINE  COUNT, 

AND  ZERO  STEP  COUNT 

WRITE  (6,9900)  T,  THETA,  FLAG,  VEL , DELTA,  WF , WP,  WK,  XH 
NUMLI.N  = NUMLIN  ♦ 1 
STEPCT  a 0.0 

RETURN 

FORMAT  STATEMENT 

990C  FORMAT  (IX,  G14.8,  2X,  G14.8,  Al,  6(2X,  G14.8)  ) 

C 

END 
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SUBROUTINE  PAGE 


I 


P 


! 


i 


C 

C 

c 

c 

c 

c 

c 

c 

c 


c 


c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 


'THIS  SUBROUTINE  PUSHES  THE  OUTPUT  TO  A NEW  PAGE  AND  PROVIQES 
ALL  HEADINGS.  INCLUDING  CURRENT  DATE  AND  TIME.  IT  IS  CALLED 
BY  THE  SUBPROCEDURE,  RWDATA,  TO  FRINT  HE AOINGS  FOR  THE  FIRST 
PAGE  OF  THE  NEW  CASE  AND  3Y  THE  SUBPROCEDJRE.  PRINTR*  FOR 
EACH  ADDITIONAL  PAGE  OF  INPUT  THAT  FOLLOWS  FOR  THAT  CASE. 

If  CALLS  NO  SUBPROCEDURES. 


INTEGER  BADXFG . 
COMMON  /FLAGS  / 

1 

COMMON  /PRINTS/ 

1 


BPFLAG,  DON  EFG . EOFLAG 

BADXFG.  BPFLAG,  DONEFG,  EOFLAG,  IERRFG,  KOUNT, 
LOADFG,  MECHFG.  M1STFG,  NCONFG 

FLAG,  MAXLIN,  N AIDA , NF,  NUMLIN,  NUMPAG,  NWAVEF, 
STEPCT , TIMNOW,  TITLE(15>,  TODAY,  TYPE 


DIMENSION  T YPS  UP ( 2 ) , T YPWAV  ( 2 » 2)  , TYPSL3(4,3),  TYPLOA(2) 

DATA  TYPSUP  / 7 H SIMPLY,  7H CL AM FEU/ , 

1 TYPWA V / :6HGENERA,  6HL  TIME,  6H  SQUAR,  6 HE  WAVE/, 

2 TYPSL3  / 8HH0RIZ CNT , 8HAL  SLAB,,  8H  EXPLCSI,  8HVE  ABOVE, 

3 BH  , 8 HVERTIC AL , 8H  WALL  , 8H  , 

4 8HH0RIZCNT , 8HAL  SLAB,,  8H  EXPLOSI,  8HVE  BELOW/, 

5 TYPLOA  / 6HLINEAR,  6H  BLAST/ 


IF  WE  ARE  COMPLETING  THE  FIRST  PAGE  FOR  CASE 
IF  ( NJMPAG.N  E. 1)  GO  TO  100 

THEN  PRINT  SUPPORT,  WAVE,  WEIGHT  TYPES  AT  BOTTOM  OF  PAGE 

WRITE  (6,9930)  TYPSUP(NF),  (T YP WAV C I,  NWAVEF ) , 1 = 1,2), 

1 (TYPSLBd,  NAIDA)  , 1 = 1,4),  TYPLO  A (LOADFG ) 

ELSE  CONTINUE 

END  IF  (END  OF  FIRST  PAGE) 

lr  WE  ARE  COMPLETING  THE  SECOND  OR  SUCCEEDING  PAGE 

100  IF  ( NJMPAG.LT. 2)  GO  TO  500 

THEN  PRINT  * MESSAGE  AT  BOTTOM  OF  PAGE 

WRITE  (6,9940) 

ELSE  CONTINUE 
END  IF  (END  OF  2ND  PAGE) 

ADD  ONE  TO  PAGE  NUMBER  AND 

PRINT  TOP  TWO  LINES  ON  NEXT  PAGE 

50(  NUMFAG  = NUMPAG  ♦ 1 

WRITE  (6,9900)  TYPE,  TODAY,  TIMNGW,  TITLE,  NUMPAG 

I?  INPUT  (FIRST  PAGE  OF  OUTPUT)  HAS  BEEN  PRINTED  FOR  THIS  CASE 
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IF  ( NJMPAG.EQ.l)  GO  TO  1000 


C 

C 


THEN  PRINT  REGULAR  OUTPUT  HEADING  LINES 
WRITE  (6,9930)  TYPSUP(NF),  ( T YP WA V ( I, NWA VEF ) , 1=1, Z) , 

(TYPSLB (I,NAIOA) , 1=1,4),  TYPLOA(LOADFG) 

WRITE  (6,9910) 

GO  TO  2000 


C 

C 


100C 


ELSE  PRINT  INPUT  HEADING  LINE 
WRITE  (6,9920) 


C 

C 

C 

C 

C 


EM D IF  (BEG  OF  1ST  PAGE) 
RESET  LINE  COUNT 


2000 


NUMLIN  = 12 
RETURN 


C 

C 

C 

C 


FORMAT  STATEMENTS 


9900  FORMAT 


9910  FORMAT 


(1H1/1H0/2 
A10,  4X» 

( 1H0/ 


1 

2 

3 

4 

5 

6 

7 

8 
9 
A 
B 
C 
D 


9920 

9930 


FORMAT 

FORMAT 


9940  FORMAT 


5X, 
1CX, 
10X, 
4X  , 
4X  , 

3 X , 
3X  , 
3X, 

3 X , 
5X, 
9X, 
6X  » 
IX, 
7X, 

(1H0/1H0  , 
(1H0,  10X, 
10X,  4 A 8 , 
( 1H0,  20  X , 
♦ELEMENT 


2X,  ♦CALCULATIONS 
A9/7X,  15A5,  43  X, 
♦TIME*, 

♦THETA*, 

♦MIDPT.  VEL.*, 

♦MIDPT.  DELTA*, 
♦PRESSURE  WORK*, 
♦PLASTIC  WORK*, 

♦KINETIC  ENERGY*, 

♦HINGE  LOCATION*/ 

* (SECONDS)*, 
♦(RADIANS)*, 

* (IN. /SEC.)*, 

♦(INCHES) *, 

3 ( 6X , * ( IN. -LBS. ) * 
♦(INCHES)*  / ) 

2X,  * IN  PUT  VALUES*  ) 

A 7,  *-SUFPORTEC* , 

10X,  A6,  * LOAD*) 

♦AN  ASTERISK  INDICATES 
HAS  FRACTURED*) 


ON  A CONCRETE 
♦PAGE*,  13) 


*,  A5,  50X, 


) 


10X,  2A6,  * FUNCTION*, 


THAT  A REINFORCING  *, 


END 
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